Add note to avoid constraints=all-bonds

grompp will now issue a note when using constraints=all-bonds with
Amber, CHARMM and OPLS-aa when using a time step smaller than 2.6 fs
and no virtual sites.

Change-Id: Ibe7a9989f202ad6072af156bd0f146a26d775c2c

diff --git a/docs/release-notes/miscellaneous.rst b/docs/release-notes/miscellaneous.rst
index aabf27d8999e87d931ac296877604fd7f21f54c8..eb5c082061aafb529af2b9d1ec13349a2e6d8ef2 100644 (file)
--- a/docs/release-notes/miscellaneous.rst
+++ b/docs/release-notes/miscellaneous.rst
@@ -1,3 +1,12 @@
 Miscellaneous
 ^^^^^^^^^^^^^
 
+grompp discourages use of constraints=all-bonds
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+Common force fields, including AMBER, CHARMM and OPLS/aa, are parametrized
+with bonds involving hydrogen constrained. Constraining all bonds should
+be avoided, for correctness. grompp now issues a note when
+constraints=all-bonds is used with these force fields when time steps
+are smaller than 2.6 fs and hydrogens are not replaced by virtual sites.
+Using constraints=h-bonds will also improve performance.
+

diff --git a/src/gromacs/gmxpreprocess/gmxcpp.cpp b/src/gromacs/gmxpreprocess/gmxcpp.cpp
index 7c0495febce7ec9552873462cebc17343b0ea9e1..2c27664586e0a1d0098ff2a8cbea84bb5438ea0f 100644 (file)
--- a/src/gromacs/gmxpreprocess/gmxcpp.cpp
+++ b/src/gromacs/gmxpreprocess/gmxcpp.cpp
@@ -684,6 +684,20 @@ int cpp_close_file(gmx_cpp_t *handlep)
     return eCPP_OK;
 }
 
+const std::string *cpp_find_define(const gmx_cpp_t   *handlep,
+                                   const std::string &defineName)
+{
+    for (const t_define &define : *(*handlep)->defines)
+    {
+        if (define.name == defineName)
+        {
+            return &define.def;
+        }
+    }
+
+    return nullptr;
+}
+
 void cpp_done(gmx_cpp_t handle)
 {
     int status = cpp_close_file(&handle);

diff --git a/src/gromacs/gmxpreprocess/gmxcpp.h b/src/gromacs/gmxpreprocess/gmxcpp.h
index ebd58bbfd5318297acad7b4f3eb40d677f1ea872..eca07f155ddb7b4cacddd864f68a0e088b063f53 100644 (file)
--- a/src/gromacs/gmxpreprocess/gmxcpp.h
+++ b/src/gromacs/gmxpreprocess/gmxcpp.h
@@ -38,6 +38,8 @@
 #ifndef GMX_GMXPREPROCESS_GMXCPP_H
 #define GMX_GMXPREPROCESS_GMXCPP_H
 
+#include <string>
+
 typedef struct gmx_cpp *gmx_cpp_t;
 
 /* The possible return codes for these functions */
@@ -70,6 +72,11 @@ int cpp_cur_linenr(const gmx_cpp_t *handlep);
  */
 int cpp_close_file(gmx_cpp_t *handlep);
 
+/* Return a pointer to the value of defineName, when present, nullptr othwerwise.
+ */
+const std::string *cpp_find_define(const gmx_cpp_t   *handlep,
+                                   const std::string &defineName);
+
 /* Clean up normal and file static data structures
  */
 void cpp_done(gmx_cpp_t handle);

diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp
index bef728f7ac844880183b5eec383d17d02826cd6b..e5bd75c893a954d1bcf7cdcbd536d1eb23b3ca02 100644 (file)
--- a/src/gromacs/gmxpreprocess/grompp.cpp
+++ b/src/gromacs/gmxpreprocess/grompp.cpp
@@ -509,6 +509,7 @@ new_status(const char *topfile, const char *topppfile, const char *confin,
     t_molinfo                  *molinfo = nullptr;
     std::vector<gmx_molblock_t> molblock;
     int                         i, nrmols, nmismatch;
+    bool                        ffParametrizedWithHBondConstraints;
     char                        buf[STRLEN];
     char                        warn_buf[STRLEN];
 
@@ -518,6 +519,7 @@ new_status(const char *topfile, const char *topppfile, const char *confin,
                        atype, &nrmols, &molinfo, intermolecular_interactions,
                        ir,
                        &molblock,
+                       &ffParametrizedWithHBondConstraints,
                        wi);
 
     sys->molblock.clear();
@@ -576,6 +578,21 @@ new_status(const char *topfile, const char *topppfile, const char *confin,
 
     gmx_mtop_finalize(sys);
 
+    /* Discourage using the, previously common, setup of applying constraints
+     * to all bonds with force fields that have been parametrized with H-bond
+     * constraints only. Do not print note with large timesteps or vsites.
+     */
+    if (opts->nshake == eshALLBONDS &&
+        ffParametrizedWithHBondConstraints &&
+        ir->delta_t < 0.0026 &&
+        gmx_mtop_ftype_count(sys, F_VSITE3FD) == 0)
+    {
+        set_warning_line(wi, "unknown", -1);
+        warning_note(wi, "You are using constraints on all bonds, whereas the forcefield "
+                     "has been parametrized only with constraints involving hydrogen atoms. "
+                     "We suggest using constraints = h-bonds instead, this will also improve performance.");
+    }
+
     /* COORDINATE file processing */
     if (bVerbose)
     {

diff --git a/src/gromacs/gmxpreprocess/topio.cpp b/src/gromacs/gmxpreprocess/topio.cpp
index a43bb17147f60f2c311d5a1e5af6fc6e277cc16a..e753c04d7d060d0333bfc74bd404452998522ea4 100644 (file)
--- a/src/gromacs/gmxpreprocess/topio.cpp
+++ b/src/gromacs/gmxpreprocess/topio.cpp
@@ -424,6 +424,7 @@ static char **read_topol(const char *infile, const char *outfile,
                          t_gromppopts *opts,
                          real        *fudgeQQ,
                          std::vector<gmx_molblock_t> *molblock,
+                         bool       *ffParametrizedWithHBondConstraints,
                          bool        bFEP,
                          bool        bZero,
                          bool        usingFullRangeElectrostatics,
@@ -903,12 +904,34 @@ static char **read_topol(const char *infile, const char *outfile,
     }
     while (!done);
     sfree(cpp_opts_return);
-    cpp_done(handle);
+
     if (out)
     {
         gmx_fio_fclose(out);
     }
 
+    /* List of GROMACS define names for force fields that have been
+     * parametrized using constraints involving hydrogens only.
+     *
+     * We should avoid hardcoded names, but this is hopefully only
+     * needed temparorily for discouraging use of constraints=all-bonds.
+     */
+    const std::array<std::string, 3> ffDefines = {
+        "_FF_AMBER",
+        "_FF_CHARMM",
+        "_FF_OPLSAA"
+    };
+    *ffParametrizedWithHBondConstraints = false;
+    for (const std::string &ffDefine : ffDefines)
+    {
+        if (cpp_find_define(&handle, ffDefine))
+        {
+            *ffParametrizedWithHBondConstraints = true;
+        }
+    }
+
+    cpp_done(handle);
+
     if (opts->couple_moltype)
     {
         if (nmol_couple == 0)
@@ -977,6 +1000,7 @@ char **do_top(bool                          bVerbose,
               t_molinfo                   **intermolecular_interactions,
               const t_inputrec             *ir,
               std::vector<gmx_molblock_t>  *molblock,
+              bool                         *ffParametrizedWithHBondConstraints,
               warninp_t                     wi)
 {
     /* Tmpfile might contain a long path */
@@ -1001,6 +1025,7 @@ char **do_top(bool                          bVerbose,
                        nrmols, molinfo, intermolecular_interactions,
                        plist, combination_rule, repulsion_power,
                        opts, fudgeQQ, molblock,
+                       ffParametrizedWithHBondConstraints,
                        ir->efep != efepNO, bZero,
                        EEL_FULL(ir->coulombtype), wi);
 

diff --git a/src/gromacs/gmxpreprocess/topio.h b/src/gromacs/gmxpreprocess/topio.h
index 2e99405b6fc8e391b78e49ec32092420bf716eb0..8d8ce585755dfa64ea3ff6a7a04d9760bc6862c8 100644 (file)
--- a/src/gromacs/gmxpreprocess/topio.h
+++ b/src/gromacs/gmxpreprocess/topio.h
@@ -69,6 +69,7 @@ char **do_top(bool                          bVerbose,
               t_molinfo                   **intermolecular_interactions,
               const t_inputrec             *ir,
               std::vector<gmx_molblock_t>  *molblock,
+              bool                         *ffParametrizedWithHBondConstraints,
               warninp_t                     wi);
 
 /* This routine expects sys->molt[m].ilist to be of size F_NRE and ordered. */