From 58a8ae3dc4c4b47ef5a5e42f9ffc1eb24ffbb17f Mon Sep 17 00:00:00 2001
From: Berk Hess <hess@kth.se>
Date: Fri, 25 Sep 2020 10:02:21 +0200
Subject: [PATCH] Fix grompp group handling with frozen system

When the whole system was frozen, grompp would attempt to write
a vector out of range.

Fixes #3683
---
 docs/release-notes/2020/2020.4.rst   |  8 ++++++++
 src/gromacs/gmxpreprocess/readir.cpp |  7 ++++---
 src/gromacs/topology/topology.h      | 11 ++++++-----
 3 files changed, 18 insertions(+), 8 deletions(-)

diff --git a/docs/release-notes/2020/2020.4.rst b/docs/release-notes/2020/2020.4.rst
index ec807dbfb7..603b9ff506 100644
--- a/docs/release-notes/2020/2020.4.rst
+++ b/docs/release-notes/2020/2020.4.rst
@@ -70,6 +70,14 @@ Only check for RDTSCP on x86 platforms
 Miscellaneous
 ^^^^^^^^^^^^^
 
+Fix crash of grompp when the whole system is frozen
+"""""""""""""""""""""""""""""""""""""""""""""""""""
+
+When the whole system would be frozen, grompp would crash with
+a segmentation fault.
+
+:issue:`3683`
+
 Fixes the unexpected change in molecule indexing in output after simulation
 """""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
 
diff --git a/src/gromacs/gmxpreprocess/readir.cpp b/src/gromacs/gmxpreprocess/readir.cpp
index 5b662c6e31..9c6d01670a 100644
--- a/src/gromacs/gmxpreprocess/readir.cpp
+++ b/src/gromacs/gmxpreprocess/readir.cpp
@@ -3268,11 +3268,12 @@ static void checkAndUpdateVcmFreezeGroupConsistency(SimulationGroups* groups,
             if (numFrozenDims == DIM)
             {
                 /* Do not remove COM motion for this fully frozen atom */
-                if (groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].empty())
+                if (groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].empty())
                 {
-                    groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval].resize(numAtoms, 0);
+                    groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval].resize(
+                            numAtoms, 0);
                 }
-                groups->groups[SimulationAtomGroupType::MassCenterVelocityRemoval][a] = vcmRestGroup;
+                groups->groupNumbers[SimulationAtomGroupType::MassCenterVelocityRemoval][a] = vcmRestGroup;
                 numFullyFrozenVcmAtoms++;
             }
             else if (numFrozenDims > 0)
diff --git a/src/gromacs/topology/topology.h b/src/gromacs/topology/topology.h
index 710bbbd2d8..1d8f9abc21 100644
--- a/src/gromacs/topology/topology.h
+++ b/src/gromacs/topology/topology.h
@@ -3,7 +3,7 @@
  *
  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  * Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2011,2014,2015,2016,2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2011,2014,2015,2016,2018,2019,2020, 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.
@@ -116,17 +116,18 @@ struct MoleculeBlockIndices
  */
 struct SimulationGroups
 {
-    //! Groups of particles
+    // TODO: collect groups and groupNumbers in a struct for each group type
+    //! Group numbers for each of the different SimulationAtomGroupType groups.
     gmx::EnumerationArray<SimulationAtomGroupType, AtomGroupIndices> groups;
     //! Names of groups, stored as pointer to the entries in the symbol table.
     std::vector<char**> groupNames;
-    //! Group numbers for the different SimulationAtomGroupType groups.
+    //! Indices into groups for each atom for each of the different SimulationAtomGroupType groups.
     gmx::EnumerationArray<SimulationAtomGroupType, std::vector<unsigned char>> groupNumbers;
 
     /*! \brief
-     * Number of group numbers for a single SimulationGroup.
+     * Number of atoms for which group numbers are stored for a single SimulationGroup.
      *
-     * \param[in] group Integer value for the group type.
+     * \param[in] group  The group type.
      */
     int numberOfGroupNumbers(SimulationAtomGroupType group) const
     {
-- 
2.22.0