Prohibit SETTLE interactions on atoms with perturbed masses

This has never been implemented, but in 2018.x and older it sometimes
managed not to crash. Now both grompp and mdrun prevent the simulation
from running and suggest an alternative.

Fixes #3959

diff --git a/docs/release-notes/2021/2021.2.rst b/docs/release-notes/2021/2021.2.rst
index 397ed3df3f41a6f7036c270851f1e532ca4a717f..fbcdafdde05b763343dca1842fc9d36ad10fdf3a 100644 (file)
--- a/docs/release-notes/2021/2021.2.rst
+++ b/docs/release-notes/2021/2021.2.rst
@@ -17,13 +17,22 @@ Fixes where mdrun could behave incorrectly
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 Removed a potential race condition with GPU update
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
 Fixed possible (but so far unobserved) race condition in coordinate copy when
 using GPU update with dipole moment calculation.
 
 :issue:`4024`
 
+Prohibited SETTLE interactions for atoms with perturbed masses
+""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+Older implementations produced varying degrees of wrong results because
+this has never been implemented. Now both ``mdrun`` and ``grompp``
+refuse to handle such a system, suggesting using normal constraints.
+
+:issue:`3959`
+
 Fixes for ``gmx`` tools
 ^^^^^^^^^^^^^^^^^^^^^^^
 

diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp
index 46a48b93e434b724e0038fa7fe0ba875430a501a..9c4fc7f255b27b597c5884446540cabf8e4f5251 100644 (file)
--- a/src/gromacs/gmxpreprocess/grompp.cpp
+++ b/src/gromacs/gmxpreprocess/grompp.cpp
@@ -4,7 +4,7 @@
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, 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.
@@ -2081,6 +2081,13 @@ int gmx_grompp(int argc, char* argv[])
     /* check masses */
     check_mol(&sys, wi);
 
+    if (haveFepPerturbedMassesInSettles(sys))
+    {
+        warning_error(wi,
+                      "SETTLE is not implemented for atoms whose mass is perturbed. "
+                      "You might instead use normal constraints.");
+    }
+
     checkForUnboundAtoms(&sys, bVerbose, wi, logger);
 
     if (EI_DYNAMICS(ir->eI) && ir->eI != eiBD)

diff --git a/src/gromacs/mdlib/constr.cpp b/src/gromacs/mdlib/constr.cpp
index 03d4783380b8566995027e32a774821549138270..32a12e0c88d0b2db90266d1361aa36a0c15631d6 100644 (file)
--- a/src/gromacs/mdlib/constr.cpp
+++ b/src/gromacs/mdlib/constr.cpp
@@ -1102,6 +1102,15 @@ Constraints::Impl::Impl(const gmx_mtop_t&     mtop_p,
 
         settled = std::make_unique<SettleData>(mtop);
 
+        // SETTLE with perturbed masses is not implemented. grompp now checks
+        // for this, but old .tpr files that did this might still exist.
+        if (haveFepPerturbedMassesInSettles(mtop))
+        {
+            gmx_fatal(FARGS,
+                      "SETTLE is not implemented for atoms whose mass is perturbed. "
+                      "You might\ninstead use normal constraints.");
+        }
+
         /* Make an atom to settle index for use in domain decomposition */
         for (size_t mt = 0; mt < mtop.moltype.size(); mt++)
         {

diff --git a/src/gromacs/topology/mtop_util.cpp b/src/gromacs/topology/mtop_util.cpp
index 789a35eced49199550f9f64bb52dc9bcce53028a..269cc533a98318e15ec30d2e1135913573d77416 100644 (file)
--- a/src/gromacs/topology/mtop_util.cpp
+++ b/src/gromacs/topology/mtop_util.cpp
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 2008,2009,2010, The GROMACS development team.
  * Copyright (c) 2012,2013,2014,2015,2016 The GROMACS development team.
- * Copyright (c) 2017,2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2017,2018,2019,2020,2021, 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.
@@ -1084,6 +1084,25 @@ bool haveFepPerturbedMasses(const gmx_mtop_t& mtop)
     return false;
 }
 
+bool haveFepPerturbedMassesInSettles(const gmx_mtop_t& mtop)
+{
+    for (const gmx_moltype_t& molt : mtop.moltype)
+    {
+        if (!molt.ilist[F_SETTLE].empty())
+        {
+            for (int a = 0; a < molt.atoms.nr; a++)
+            {
+                const t_atom& atom = molt.atoms.atom[a];
+                if (atom.m != atom.mB)
+                {
+                    return true;
+                }
+            }
+        }
+    }
+    return false;
+}
+
 bool havePerturbedConstraints(const gmx_mtop_t& mtop)
 {
     // This code assumes that all perturbed constraints parameters are actually used

diff --git a/src/gromacs/topology/mtop_util.h b/src/gromacs/topology/mtop_util.h
index 7d0527e654cfa4f844f39acb30ed06111e7b1a91..13c69db4781f9e7aa254b58c6cdcbf301ace7aaa 100644 (file)
--- a/src/gromacs/topology/mtop_util.h
+++ b/src/gromacs/topology/mtop_util.h
@@ -4,7 +4,7 @@
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
  * Copyright (c) 2012,2013,2014,2015,2016 by the GROMACS development team.
- * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020,2021, 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.
@@ -283,6 +283,9 @@ bool haveFepPerturbedNBInteractions(const gmx_mtop_t& mtop);
 //! Checks whether masses are perturbed for free-energy calculations
 bool haveFepPerturbedMasses(const gmx_mtop_t& mtop);
 
+//! Checks whether masses are perturbed for free-energy calculations in SETTLE interactions
+bool haveFepPerturbedMassesInSettles(const gmx_mtop_t& mtop);
+
 //! Checks whether constraints are perturbed for free-energy calculations
 bool havePerturbedConstraints(const gmx_mtop_t& mtop);