Fix grompp constraints=all-angles with constraints

When using the mdp option constraints=all-angles, angles involving
bonds supplied as constraints in the topology would be removed,
but not replaced by angle constraints.

Fixes #3067

Change-Id: I7db390785a23c67abb30ae98e12272f88f5d5af7

diff --git a/docs/release-notes/2019/2019.4.rst b/docs/release-notes/2019/2019.4.rst
index 3d3f59af7f316d73185c903404cccf40f2214c6d..b58eb1d6a49b007d773ece9f854bb94e50916567 100644 (file)
--- a/docs/release-notes/2019/2019.4.rst
+++ b/docs/release-notes/2019/2019.4.rst
@@ -37,6 +37,15 @@ Fix segmentation fault in grompp and mdrun with cosine COM pulling
 Fixes for ``gmx`` tools
 ^^^^^^^^^^^^^^^^^^^^^^^
 
+Fix grompp not adding angle constraints between constraints
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+When using the mdp option constraints=all-angles, angles involving
+bonds supplied as constraints in the topology would be removed,
+but not replaced by angle constraints.
+
+:issue:`3067`
+
 Fix bug in gmx xpm2ps
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 

diff --git a/src/gromacs/gmxpreprocess/topshake.cpp b/src/gromacs/gmxpreprocess/topshake.cpp
index cf48da1d8a85c6ac6e23df7959759edaec85fbba..d8ee49e80d77e68c0ad13dbc89158db13fb952e0 100644 (file)
--- a/src/gromacs/gmxpreprocess/topshake.cpp
+++ b/src/gromacs/gmxpreprocess/topshake.cpp
@@ -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,2018, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2018,2019, 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.
@@ -139,7 +139,7 @@ void make_shake (t_params plist[], t_atoms *atoms, int nshake)
              */
             for (ftype = 0; (ftype < F_NRE); ftype++)
             {
-                if (interaction_function[ftype].flags & IF_BTYPE)
+                if (interaction_function[ftype].flags & IF_CHEMBOND)
                 {
                     bonds = &(plist[ftype]);
 
@@ -193,6 +193,12 @@ void make_shake (t_params plist[], t_atoms *atoms, int nshake)
                                     }
                                     if (bFound)
                                     {
+                                        if (ftype == F_CONNBONDS ||
+                                            ftype_a == F_CONNBONDS)
+                                        {
+                                            gmx_fatal(FARGS, "Can not constrain all angles when they involved bonds of type %s",
+                                                      interaction_function[F_CONNBONDS].longname);
+                                        }
                                         /* apply law of cosines */
                                         p.c0() = std::sqrt( b_ij*b_ij + b_jk*b_jk -
                                                             2.0*b_ij*b_jk*cos(DEG2RAD*ang->c0()) );