From 63b3d45c6f6f46b3c5843cae514819daea8f5ee4 Mon Sep 17 00:00:00 2001
From: Paul Bauer <paul.bauer.q@gmail.com>
Date: Thu, 9 May 2019 16:17:37 +0200
Subject: [PATCH] Fix membrane embedding

Refactoring forgot to allocate some data, causing a segmentation fault.
More refactoring used wrong data, causing a crash further on.
Even more refactoring lead to invalid memory access.

Fixes #2947

Change-Id: I61ad3125102b50c0338e0935a75cd7b1de95bc3f
---
 docs/release-notes/2019/2019.3.rst |  4 ++++
 src/gromacs/mdlib/membed.cpp       | 22 ++++++++++++++--------
 2 files changed, 18 insertions(+), 8 deletions(-)

diff --git a/docs/release-notes/2019/2019.3.rst b/docs/release-notes/2019/2019.3.rst
index bfea19a864..ad12d19dcd 100644
--- a/docs/release-notes/2019/2019.3.rst
+++ b/docs/release-notes/2019/2019.3.rst
@@ -24,6 +24,10 @@ a crash due to a divide by zero error. Fixed by introducing a check.
 
 :issue:`2917`
 
+Fix segmentation fault when using membrane embedding
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+:issue:`2947`
 
 Fixes for ``gmx`` tools
 ^^^^^^^^^^^^^^^^^^^^^^^
diff --git a/src/gromacs/mdlib/membed.cpp b/src/gromacs/mdlib/membed.cpp
index c564756966..ad70bf90ad 100644
--- a/src/gromacs/mdlib/membed.cpp
+++ b/src/gromacs/mdlib/membed.cpp
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2010,2011,2012,2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
+ * Copyright (c) 2010,2011,2012,2013,2014,2015,2016,2017,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.
@@ -530,6 +530,7 @@ static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc
     r_min_rad = probe_rad*probe_rad;
     gmx::RangePartitioning molecules = gmx_mtop_molecules(*mtop);
     snew(rm_p->block, molecules.numBlocks());
+    snew(rm_p->mol, molecules.numBlocks());
     nrm    = nupper = 0;
     nlower = low_up_rm;
     for (i = 0; i < ins_at->nr; i++)
@@ -689,7 +690,7 @@ static void rm_group(gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state
     for (int i = 0; i < rm_p->nr; i++)
     {
         mol_id = rm_p->mol[i];
-        at     = molecules.block(mol_id).size();
+        at     = molecules.block(mol_id).begin();
         block  = rm_p->block[i];
         mtop->molblock[block].nmol--;
         for (j = 0; j < mtop->moltype[mtop->molblock[block].type].atoms.nr; j++)
@@ -700,23 +701,27 @@ static void rm_group(gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state
     }
 
     mtop->natoms    -= n;
-    state_change_natoms(state, state->natoms - n);
-    snew(x_tmp, state->natoms);
-    snew(v_tmp, state->natoms);
+    /* We cannot change the size of the state datastructures here
+     * because we still access the coordinate arrays for all positions
+     * before removing the molecules we want to remove.
+     */
+    const int newStateAtomNumber = state->natoms - n;
+    snew(x_tmp, newStateAtomNumber);
+    snew(v_tmp, newStateAtomNumber);
 
     for (int i = 0; i < egcNR; i++)
     {
         if (groups->grpnr[i] != nullptr)
         {
-            groups->ngrpnr[i] = state->natoms;
-            snew(new_egrp[i], state->natoms);
+            groups->ngrpnr[i] = newStateAtomNumber;
+            snew(new_egrp[i], newStateAtomNumber);
         }
     }
 
     auto x = makeArrayRef(state->x);
     auto v = makeArrayRef(state->v);
     rm = 0;
-    for (int i = 0; i < state->natoms+n; i++)
+    for (int i = 0; i < state->natoms; i++)
     {
         bRM = FALSE;
         for (j = 0; j < n; j++)
@@ -759,6 +764,7 @@ static void rm_group(gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state
             }
         }
     }
+    state_change_natoms(state, newStateAtomNumber);
     for (int i = 0; i < state->natoms; i++)
     {
         copy_rvec(x_tmp[i], x[i]);
-- 
2.11.4.GIT