From 18840da5f2766cf4686c139db5138bdddcb74c47 Mon Sep 17 00:00:00 2001
From: Carlo Camilloni <carlo.camilloni@gmail.com>
Date: Sun, 7 Aug 2016 11:35:13 +0200
Subject: [PATCH] Gromacs 2016 patch

---
 .../src/gromacs/CMakeLists.txt                |   5 +
 .../src/gromacs/CMakeLists.txt.preplumed      |   5 +
 .../src/gromacs/mdlib/force.cpp               |   0
 .../src/gromacs/mdlib/force.cpp.preplumed     |   0
 .../src/gromacs/mdlib/minimize.cpp            |   5 +
 .../src/gromacs/mdlib/minimize.cpp.preplumed  |   5 +
 .../src/programs/mdrun/md.cpp                 | 203 +++++++++---------
 .../src/programs/mdrun/md.cpp.preplumed       | 203 +++++++++---------
 .../src/programs/mdrun/mdrun.cpp              |  24 ++-
 .../src/programs/mdrun/mdrun.cpp.preplumed    |  24 ++-
 .../src/programs/mdrun/repl_ex.cpp            |   0
 .../src/programs/mdrun/repl_ex.cpp.preplumed  |   0
 .../src/programs/mdrun/runner.cpp             |  25 ++-
 .../src/programs/mdrun/runner.cpp.preplumed   |  25 ++-
 14 files changed, 314 insertions(+), 210 deletions(-)
 rename patches/{gromacs-2016-beta2.diff => gromacs-2016.diff}/src/gromacs/CMakeLists.txt (96%)
 rename patches/{gromacs-2016-beta2.diff => gromacs-2016.diff}/src/gromacs/CMakeLists.txt.preplumed (96%)
 rename patches/{gromacs-2016-beta2.diff => gromacs-2016.diff}/src/gromacs/mdlib/force.cpp (100%)
 rename patches/{gromacs-2016-beta2.diff => gromacs-2016.diff}/src/gromacs/mdlib/force.cpp.preplumed (100%)
 rename patches/{gromacs-2016-beta2.diff => gromacs-2016.diff}/src/gromacs/mdlib/minimize.cpp (99%)
 rename patches/{gromacs-2016-beta2.diff => gromacs-2016.diff}/src/gromacs/mdlib/minimize.cpp.preplumed (99%)
 rename patches/{gromacs-2016-beta2.diff => gromacs-2016.diff}/src/programs/mdrun/md.cpp (92%)
 rename patches/{gromacs-2016-beta2.diff => gromacs-2016.diff}/src/programs/mdrun/md.cpp.preplumed (92%)
 rename patches/{gromacs-2016-beta2.diff => gromacs-2016.diff}/src/programs/mdrun/mdrun.cpp (97%)
 rename patches/{gromacs-2016-beta2.diff => gromacs-2016.diff}/src/programs/mdrun/mdrun.cpp.preplumed (96%)
 rename patches/{gromacs-2016-beta2.diff => gromacs-2016.diff}/src/programs/mdrun/repl_ex.cpp (100%)
 rename patches/{gromacs-2016-beta2.diff => gromacs-2016.diff}/src/programs/mdrun/repl_ex.cpp.preplumed (100%)
 rename patches/{gromacs-2016-beta2.diff => gromacs-2016.diff}/src/programs/mdrun/runner.cpp (98%)
 rename patches/{gromacs-2016-beta2.diff => gromacs-2016.diff}/src/programs/mdrun/runner.cpp.preplumed (98%)

diff --git a/patches/gromacs-2016-beta2.diff/src/gromacs/CMakeLists.txt b/patches/gromacs-2016.diff/src/gromacs/CMakeLists.txt
similarity index 96%
rename from patches/gromacs-2016-beta2.diff/src/gromacs/CMakeLists.txt
rename to patches/gromacs-2016.diff/src/gromacs/CMakeLists.txt
index a9f3b4fd3..dfb7260df 100644
--- a/patches/gromacs-2016-beta2.diff/src/gromacs/CMakeLists.txt
+++ b/patches/gromacs-2016.diff/src/gromacs/CMakeLists.txt
@@ -126,6 +126,10 @@ list(APPEND LIBGROMACS_SOURCES ${THREAD_MPI_SOURCES})
 list(APPEND LIBGROMACS_SOURCES ${TNG_SOURCES})
 tng_set_source_properties(WITH_ZLIB ${HAVE_ZLIB})
 
+get_lmfit_properties(LMFIT_SOURCES LMFIT_LIBRARIES_TO_LINK LMFIT_INCLUDE_DIRECTORY LMFIT_INCLUDE_DIR_ORDER)
+include_directories(${LMFIT_INCLUDE_DIR_ORDER} SYSTEM "${LMFIT_INCLUDE_DIRECTORY}")
+list(APPEND LIBGROMACS_SOURCES ${LMFIT_SOURCES})
+
 configure_file(version.h.cmakein version.h)
 gmx_install_headers(
     analysisdata.h
@@ -183,6 +187,7 @@ target_link_libraries(libgromacs
                       ${TNG_IO_LIBRARIES}
                       ${FFT_LIBRARIES} ${LINEAR_ALGEBRA_LIBRARIES}
                       ${XML_LIBRARIES}
+                      ${LMFIT_LIBRARIES_TO_LINK}
                       ${THREAD_LIB} ${GMX_SHARED_LINKER_FLAGS} ${OPENCL_LIBRARIES}
                       ${GMX_STDLIB_LIBRARIES})
 set_target_properties(libgromacs PROPERTIES
diff --git a/patches/gromacs-2016-beta2.diff/src/gromacs/CMakeLists.txt.preplumed b/patches/gromacs-2016.diff/src/gromacs/CMakeLists.txt.preplumed
similarity index 96%
rename from patches/gromacs-2016-beta2.diff/src/gromacs/CMakeLists.txt.preplumed
rename to patches/gromacs-2016.diff/src/gromacs/CMakeLists.txt.preplumed
index ab9d64df4..fd321f1e1 100644
--- a/patches/gromacs-2016-beta2.diff/src/gromacs/CMakeLists.txt.preplumed
+++ b/patches/gromacs-2016.diff/src/gromacs/CMakeLists.txt.preplumed
@@ -124,6 +124,10 @@ list(APPEND LIBGROMACS_SOURCES ${THREAD_MPI_SOURCES})
 list(APPEND LIBGROMACS_SOURCES ${TNG_SOURCES})
 tng_set_source_properties(WITH_ZLIB ${HAVE_ZLIB})
 
+get_lmfit_properties(LMFIT_SOURCES LMFIT_LIBRARIES_TO_LINK LMFIT_INCLUDE_DIRECTORY LMFIT_INCLUDE_DIR_ORDER)
+include_directories(${LMFIT_INCLUDE_DIR_ORDER} SYSTEM "${LMFIT_INCLUDE_DIRECTORY}")
+list(APPEND LIBGROMACS_SOURCES ${LMFIT_SOURCES})
+
 configure_file(version.h.cmakein version.h)
 gmx_install_headers(
     analysisdata.h
@@ -179,6 +183,7 @@ target_link_libraries(libgromacs
                       ${TNG_IO_LIBRARIES}
                       ${FFT_LIBRARIES} ${LINEAR_ALGEBRA_LIBRARIES}
                       ${XML_LIBRARIES}
+                      ${LMFIT_LIBRARIES_TO_LINK}
                       ${THREAD_LIB} ${GMX_SHARED_LINKER_FLAGS} ${OPENCL_LIBRARIES}
                       ${GMX_STDLIB_LIBRARIES})
 set_target_properties(libgromacs PROPERTIES
diff --git a/patches/gromacs-2016-beta2.diff/src/gromacs/mdlib/force.cpp b/patches/gromacs-2016.diff/src/gromacs/mdlib/force.cpp
similarity index 100%
rename from patches/gromacs-2016-beta2.diff/src/gromacs/mdlib/force.cpp
rename to patches/gromacs-2016.diff/src/gromacs/mdlib/force.cpp
diff --git a/patches/gromacs-2016-beta2.diff/src/gromacs/mdlib/force.cpp.preplumed b/patches/gromacs-2016.diff/src/gromacs/mdlib/force.cpp.preplumed
similarity index 100%
rename from patches/gromacs-2016-beta2.diff/src/gromacs/mdlib/force.cpp.preplumed
rename to patches/gromacs-2016.diff/src/gromacs/mdlib/force.cpp.preplumed
diff --git a/patches/gromacs-2016-beta2.diff/src/gromacs/mdlib/minimize.cpp b/patches/gromacs-2016.diff/src/gromacs/mdlib/minimize.cpp
similarity index 99%
rename from patches/gromacs-2016-beta2.diff/src/gromacs/mdlib/minimize.cpp
rename to patches/gromacs-2016.diff/src/gromacs/mdlib/minimize.cpp
index b60b80991..6282d144c 100644
--- a/patches/gromacs-2016-beta2.diff/src/gromacs/mdlib/minimize.cpp
+++ b/patches/gromacs-2016.diff/src/gromacs/mdlib/minimize.cpp
@@ -1089,6 +1089,7 @@ namespace gmx
                            gmx_edsam_t ed,
                            t_forcerec *fr,
                            int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
+                           gmx_membed_t gmx_unused *membed,
                            real cpt_period, real max_hours,
                            int imdport,
                            unsigned long Flags,
@@ -1108,6 +1109,7 @@ double do_cg(FILE *fplog, t_commrec *cr,
              gmx_edsam_t gmx_unused ed,
              t_forcerec *fr,
              int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed,
+             gmx_membed_t gmx_unused *membed,
              real gmx_unused cpt_period, real gmx_unused max_hours,
              int imdport,
              unsigned long gmx_unused Flags,
@@ -1754,6 +1756,7 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
                 gmx_edsam_t gmx_unused ed,
                 t_forcerec *fr,
                 int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed,
+                gmx_membed_t gmx_unused *membed,
                 real gmx_unused cpt_period, real gmx_unused max_hours,
                 int imdport,
                 unsigned long gmx_unused Flags,
@@ -2581,6 +2584,7 @@ double do_steep(FILE *fplog, t_commrec *cr,
                 gmx_edsam_t gmx_unused  ed,
                 t_forcerec *fr,
                 int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed,
+                gmx_membed_t gmx_unused *membed,
                 real gmx_unused cpt_period, real gmx_unused max_hours,
                 int imdport,
                 unsigned long gmx_unused Flags,
@@ -2848,6 +2852,7 @@ double do_nm(FILE *fplog, t_commrec *cr,
              gmx_edsam_t  gmx_unused ed,
              t_forcerec *fr,
              int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed,
+             gmx_membed_t gmx_unused *membed,
              real gmx_unused cpt_period, real gmx_unused max_hours,
              int imdport,
              unsigned long gmx_unused Flags,
diff --git a/patches/gromacs-2016-beta2.diff/src/gromacs/mdlib/minimize.cpp.preplumed b/patches/gromacs-2016.diff/src/gromacs/mdlib/minimize.cpp.preplumed
similarity index 99%
rename from patches/gromacs-2016-beta2.diff/src/gromacs/mdlib/minimize.cpp.preplumed
rename to patches/gromacs-2016.diff/src/gromacs/mdlib/minimize.cpp.preplumed
index 55200ce0a..3890e24e3 100644
--- a/patches/gromacs-2016-beta2.diff/src/gromacs/mdlib/minimize.cpp.preplumed
+++ b/patches/gromacs-2016.diff/src/gromacs/mdlib/minimize.cpp.preplumed
@@ -1015,6 +1015,7 @@ namespace gmx
                            gmx_edsam_t ed,
                            t_forcerec *fr,
                            int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
+                           gmx_membed_t gmx_unused *membed,
                            real cpt_period, real max_hours,
                            int imdport,
                            unsigned long Flags,
@@ -1034,6 +1035,7 @@ double do_cg(FILE *fplog, t_commrec *cr,
              gmx_edsam_t gmx_unused ed,
              t_forcerec *fr,
              int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed,
+             gmx_membed_t gmx_unused *membed,
              real gmx_unused cpt_period, real gmx_unused max_hours,
              int imdport,
              unsigned long gmx_unused Flags,
@@ -1680,6 +1682,7 @@ double do_lbfgs(FILE *fplog, t_commrec *cr,
                 gmx_edsam_t gmx_unused ed,
                 t_forcerec *fr,
                 int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed,
+                gmx_membed_t gmx_unused *membed,
                 real gmx_unused cpt_period, real gmx_unused max_hours,
                 int imdport,
                 unsigned long gmx_unused Flags,
@@ -2507,6 +2510,7 @@ double do_steep(FILE *fplog, t_commrec *cr,
                 gmx_edsam_t gmx_unused  ed,
                 t_forcerec *fr,
                 int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed,
+                gmx_membed_t gmx_unused *membed,
                 real gmx_unused cpt_period, real gmx_unused max_hours,
                 int imdport,
                 unsigned long gmx_unused Flags,
@@ -2774,6 +2778,7 @@ double do_nm(FILE *fplog, t_commrec *cr,
              gmx_edsam_t  gmx_unused ed,
              t_forcerec *fr,
              int gmx_unused repl_ex_nst, int gmx_unused repl_ex_nex, int gmx_unused repl_ex_seed,
+             gmx_membed_t gmx_unused *membed,
              real gmx_unused cpt_period, real gmx_unused max_hours,
              int imdport,
              unsigned long gmx_unused Flags,
diff --git a/patches/gromacs-2016-beta2.diff/src/programs/mdrun/md.cpp b/patches/gromacs-2016.diff/src/programs/mdrun/md.cpp
similarity index 92%
rename from patches/gromacs-2016-beta2.diff/src/programs/mdrun/md.cpp
rename to patches/gromacs-2016.diff/src/programs/mdrun/md.cpp
index 19cea9f69..90633b424 100644
--- a/patches/gromacs-2016-beta2.diff/src/programs/mdrun/md.cpp
+++ b/patches/gromacs-2016.diff/src/programs/mdrun/md.cpp
@@ -77,13 +77,13 @@
 #include "gromacs/mdlib/mdebin.h"
 #include "gromacs/mdlib/mdoutf.h"
 #include "gromacs/mdlib/mdrun.h"
-#include "gromacs/mdlib/mdrun_signalling.h"
 #include "gromacs/mdlib/nb_verlet.h"
 #include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
 #include "gromacs/mdlib/ns.h"
 #include "gromacs/mdlib/shellfc.h"
 #include "gromacs/mdlib/sighandler.h"
 #include "gromacs/mdlib/sim_util.h"
+#include "gromacs/mdlib/simulationsignal.h"
 #include "gromacs/mdlib/tgroup.h"
 #include "gromacs/mdlib/trajectory_writing.h"
 #include "gromacs/mdlib/update.h"
@@ -131,6 +131,8 @@ extern plumed plumedmain;
 #include "corewrap.h"
 #endif
 
+using gmx::SimulationSignaller;
+
 /*! \brief Check whether bonded interactions are missing, if appropriate
  *
  * \param[in]    fplog                                  Log file pointer
@@ -224,6 +226,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                   t_nrnb *nrnb, gmx_wallcycle_t wcycle,
                   gmx_edsam_t ed, t_forcerec *fr,
                   int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
+                  gmx_membed_t *membed,
                   real cpt_period, real max_hours,
                   int imdport,
                   unsigned long Flags,
@@ -236,7 +239,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     gmx_bool        bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
     gmx_bool        bNS, bNStList, bSimAnn, bStopCM, bRerunMD,
                     bFirstStep, startingFromCheckpoint, bInitStep, bLastStep = FALSE,
-                    bBornRadii;
+                    bBornRadii, bUsingEnsembleRestraints;
     gmx_bool          bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
     gmx_bool          do_ene, do_log, do_verbose, bRerunWarnNoV = TRUE,
                       bForceUpdate = FALSE, bCPT;
@@ -259,7 +262,6 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     gmx_global_stat_t gstat;
     gmx_update_t     *upd   = NULL;
     t_graph          *graph = NULL;
-    gmx_signalling_t  gs;
     gmx_groups_t     *groups;
     gmx_ekindata_t   *ekind;
     gmx_shellfc_t    *shellfc;
@@ -280,9 +282,8 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     int             **trotter_seq;
     char              sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
     int               handled_stop_condition = gmx_stop_cond_none; /* compare to get_stop_condition*/
-    gmx_int64_t       multisim_nsteps        = -1;                 /* number of steps to do  before first multisim
-                                                                          simulation stops. If equal to zero, don't
-                                                                          communicate any more between multisims.*/
+
+
     /* PME load balancing data for GPU kernels */
     pme_load_balancing_t *pme_loadbal      = NULL;
     gmx_bool              bPMETune         = FALSE;
@@ -290,7 +291,6 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
 
     /* Interactive MD */
     gmx_bool          bIMDstep = FALSE;
-    gmx_membed_t     *membed   = NULL;
 
     /* PLUMED */
     int plumedNeedsEnergy=0;
@@ -308,8 +308,13 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
        code. So we do that alongside the first global energy reduction
        after a new DD is made. These variables handle whether the
        check happens, and the result it returns. */
-    bool shouldCheckNumberOfBondedInteractions = false;
-    int  totalNumberOfBondedInteractions       = -1;
+    bool              shouldCheckNumberOfBondedInteractions = false;
+    int               totalNumberOfBondedInteractions       = -1;
+
+    SimulationSignals signals;
+    // Most global communnication stages don't propagate mdrun
+    // signals, and will use this object to achieve that.
+    SimulationSignaller nullSignaller(nullptr, nullptr, false, false);
 
     /* Check for special mdrun options */
     bRerunMD = (Flags & MD_RERUN);
@@ -348,18 +353,6 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     }
     groups = &top_global->groups;
 
-    if (opt2bSet("-membed", nfile, fnm))
-    {
-        if (MASTER(cr))
-        {
-            fprintf(stderr, "Initializing membed");
-        }
-        /* Note that membed cannot work in parallel because mtop is
-         * changed here. Fix this if we ever want to make it run with
-         * multiple ranks. */
-        membed = init_membed(fplog, nfile, fnm, top_global, ir, state_global, cr, &cpt_period);
-    }
-
     if (ir->eSwapCoords != eswapNO)
     {
         /* Initialize ion swapping code */
@@ -605,7 +598,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     bSumEkinhOld = FALSE;
     compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
                     NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
-                    constr, NULL, FALSE, state->box,
+                    constr, &nullSignaller, state->box,
                     &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags
                     | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
     checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
@@ -621,7 +614,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
 
         compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
                         NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
-                        constr, NULL, FALSE, state->box,
+                        constr, &nullSignaller, state->box,
                         NULL, &bSumEkinhOld,
                         cglo_flags &~(CGLO_STOPCM | CGLO_PRESSURE));
     }
@@ -818,22 +811,46 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     bSumEkinhOld     = FALSE;
     bExchanged       = FALSE;
     bNeedRepartition = FALSE;
+    // TODO This implementation of ensemble orientation restraints is nasty because
+    // a user can't just do multi-sim with single-sim orientation restraints.
+    bUsingEnsembleRestraints = (fcd->disres.nsystems > 1) || (cr->ms && fcd->orires.nr);
 
-    init_global_signals(&gs, cr, ir, repl_ex_nst);
+    {
+        // Replica exchange and ensemble restraints need all
+        // simulations to remain synchronized, so they need
+        // checkpoints and stop conditions to act on the same step, so
+        // the propagation of such signals must take place between
+        // simulations, not just within simulations.
+        bool checkpointIsLocal    = (repl_ex_nst <= 0) && !bUsingEnsembleRestraints;
+        bool stopConditionIsLocal = (repl_ex_nst <= 0) && !bUsingEnsembleRestraints;
+        bool resetCountersIsLocal = true;
+        signals[eglsCHKPT]         = SimulationSignal(checkpointIsLocal);
+        signals[eglsSTOPCOND]      = SimulationSignal(stopConditionIsLocal);
+        signals[eglsRESETCOUNTERS] = SimulationSignal(resetCountersIsLocal);
+    }
 
     step     = ir->init_step;
     step_rel = 0;
 
-    if (MULTISIM(cr) && (repl_ex_nst <= 0 ))
+    // TODO extract this to new multi-simulation module
+    if (MASTER(cr) && MULTISIM(cr) && (repl_ex_nst <= 0 ))
     {
-        /* check how many steps are left in other sims */
-        multisim_nsteps = get_multisim_nsteps(cr, ir->nsteps);
+        if (!multisim_int_all_are_equal(cr->ms, ir->nsteps))
+        {
+            md_print_info(cr, fplog,
+                          "Note: The number of steps is not consistent across multi simulations,\n"
+                          "but we are proceeding anyway!\n");
+        }
+        if (!multisim_int_all_are_equal(cr->ms, ir->init_step))
+        {
+            md_print_info(cr, fplog,
+                          "Note: The initial step is not consistent across multi simulations,\n"
+                          "but we are proceeding anyway!\n");
+        }
     }
 
-
     /* and stop now if we should */
-    bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps) ||
-                 ((multisim_nsteps >= 0) && (step_rel >= multisim_nsteps )));
+    bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
     while (!bLastStep)
     {
 
@@ -876,6 +893,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             t         = t0 + step*ir->delta_t;
         }
 
+        // TODO Refactor this, so that nstfep does not need a default value of zero
         if (ir->efep != efepNO || ir->bSimTemp)
         {
             /* find and set the current lambdas.  If rerunning, we either read in a state, or a lambda value,
@@ -961,7 +979,6 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         {
             /* for rerun MD always do Neighbour Searching */
             bNS      = (bFirstStep || ir->nstlist != 0);
-            bNStList = bNS;
         }
         else
         {
@@ -969,30 +986,9 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
         }
 
-        /* check whether we should stop because another simulation has
-           stopped. */
-        if (MULTISIM(cr))
-        {
-            if ( (multisim_nsteps >= 0) &&  (step_rel >= multisim_nsteps)  &&
-                 (multisim_nsteps != ir->nsteps) )
-            {
-                if (bNS)
-                {
-                    if (MASTER(cr))
-                    {
-                        fprintf(stderr,
-                                "Stopping simulation %d because another one has finished\n",
-                                cr->ms->sim);
-                    }
-                    bLastStep         = TRUE;
-                    gs.sig[eglsCHKPT] = 1;
-                }
-            }
-        }
-
-        /* < 0 means stop after this step, > 0 means stop at next NS step */
-        if ( (gs.set[eglsSTOPCOND] < 0) ||
-             ( (gs.set[eglsSTOPCOND] > 0) && (bNStList || ir->nstlist == 0) ) )
+        /* < 0 means stop at next step, > 0 means stop at next NS step */
+        if ( (signals[eglsSTOPCOND].set < 0) ||
+             ( (signals[eglsSTOPCOND].set > 0 ) && ( bNS || ir->nstlist == 0)))
         {
             bLastStep = TRUE;
         }
@@ -1077,7 +1073,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             /* This may not be quite working correctly yet . . . . */
             compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
                             wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
-                            constr, NULL, FALSE, state->box,
+                            constr, &nullSignaller, state->box,
                             &totalNumberOfBondedInteractions, &bSumEkinhOld,
                             CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
             checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
@@ -1091,12 +1087,12 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
          * or at the last step (but not when we do not want confout),
          * but never at the first step or with rerun.
          */
-        bCPT = (((gs.set[eglsCHKPT] && (bNS || ir->nstlist == 0)) ||
+        bCPT = (((signals[eglsCHKPT].set && (bNS || ir->nstlist == 0)) ||
                  (bLastStep && (Flags & MD_CONFOUT))) &&
                 step > ir->init_step && !bRerunMD);
         if (bCPT)
         {
-            gs.set[eglsCHKPT] = 0;
+            signals[eglsCHKPT].set = 0;
         }
 
         /* Determine the energy and pressure:
@@ -1271,7 +1267,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                 wallcycle_stop(wcycle, ewcUPDATE);
                 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
-                                constr, NULL, FALSE, state->box,
+                                constr, &nullSignaller, state->box,
                                 &totalNumberOfBondedInteractions, &bSumEkinhOld,
                                 (bGStat ? CGLO_GSTAT : 0)
                                 | CGLO_ENERGY
@@ -1319,7 +1315,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                     /* This may not be quite working correctly yet . . . . */
                     compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
                                     wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
-                                    constr, NULL, FALSE, state->box,
+                                    constr, &nullSignaller, state->box,
                                     NULL, &bSumEkinhOld,
                                     CGLO_GSTAT | CGLO_TEMPERATURE);
                     wallcycle_start(wcycle, ewcUPDATE);
@@ -1398,18 +1394,18 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         {
             int nsteps_stop = -1;
 
-            /* this is just make gs.sig compatible with the hack
-               of sending signals around by MPI_Reduce with together with
+            /* this just makes signals[].sig compatible with the hack
+               of sending signals around by MPI_Reduce together with
                other floats */
             if (gmx_get_stop_condition() == gmx_stop_cond_next_ns)
             {
-                gs.sig[eglsSTOPCOND] = 1;
-                nsteps_stop          = std::max(ir->nstlist, 2*nstglobalcomm);
+                signals[eglsSTOPCOND].sig = 1;
+                nsteps_stop               = std::max(ir->nstlist, 2*nstglobalcomm);
             }
             if (gmx_get_stop_condition() == gmx_stop_cond_next)
             {
-                gs.sig[eglsSTOPCOND] = -1;
-                nsteps_stop          = nstglobalcomm + 1;
+                signals[eglsSTOPCOND].sig = -1;
+                nsteps_stop               = nstglobalcomm + 1;
             }
             if (fplog)
             {
@@ -1426,10 +1422,10 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         }
         else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
                  (max_hours > 0 && elapsed_time > max_hours*60.0*60.0*0.99) &&
-                 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
+                 signals[eglsSTOPCOND].sig == 0 && signals[eglsSTOPCOND].set == 0)
         {
             /* Signal to terminate the run */
-            gs.sig[eglsSTOPCOND] = 1;
+            signals[eglsSTOPCOND].sig = 1;
             if (fplog)
             {
                 fprintf(fplog, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step, sbuf), max_hours*0.99);
@@ -1441,7 +1437,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             elapsed_time > max_hours*60.0*60.0*0.495)
         {
             /* Set flag that will communicate the signal to all ranks in the simulation */
-            gs.sig[eglsRESETCOUNTERS] = 1;
+            signals[eglsRESETCOUNTERS].sig = 1;
         }
 
         /* In parallel we only have to check for checkpointing in steps
@@ -1452,9 +1448,9 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                            cpt_period >= 0 &&
                            (cpt_period == 0 ||
                             elapsed_time >= nchkpt*cpt_period*60.0)) &&
-            gs.set[eglsCHKPT] == 0)
+            signals[eglsCHKPT].set == 0)
         {
-            gs.sig[eglsCHKPT] = 1;
+            signals[eglsCHKPT].sig = 1;
         }
 
         /* #########   START SECOND UPDATE STEP ################# */
@@ -1544,7 +1540,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                 /* just compute the kinetic energy at the half step to perform a trotter step */
                 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
-                                constr, NULL, FALSE, lastbox,
+                                constr, &nullSignaller, lastbox,
                                 NULL, &bSumEkinhOld,
                                 (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
                                 );
@@ -1620,26 +1616,40 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
          * non-communication steps, but we need to calculate
          * the kinetic energy one step before communication.
          */
-        if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)))
         {
-            compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
-                            wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
-                            constr, &gs,
-                            (step_rel % gs.nstms == 0) &&
-                            (multisim_nsteps < 0 || (step_rel < multisim_nsteps)),
-                            lastbox,
-                            &totalNumberOfBondedInteractions, &bSumEkinhOld,
-                            (bGStat ? CGLO_GSTAT : 0)
-                            | (!EI_VV(ir->eI) || bRerunMD ? CGLO_ENERGY : 0)
-                            | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
-                            | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
-                            | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
-                            | CGLO_CONSTRAINT
-                            | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
-                            );
-            checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
-                                            top_global, top, state,
-                                            &shouldCheckNumberOfBondedInteractions);
+            // Organize to do inter-simulation signalling on steps if
+            // and when algorithms require it.
+            bool doInterSimSignal = (!bFirstStep && bDoReplEx) || bUsingEnsembleRestraints;
+
+            if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)) || doInterSimSignal)
+            {
+                // Since we're already communicating at this step, we
+                // can propagate intra-simulation signals. Note that
+                // check_nstglobalcomm has the responsibility for
+                // choosing the value of nstglobalcomm that is one way
+                // bGStat becomes true, so we can't get into a
+                // situation where e.g. checkpointing can't be
+                // signalled.
+                bool                doIntraSimSignal = true;
+                SimulationSignaller signaller(&signals, cr, doInterSimSignal, doIntraSimSignal);
+
+                compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
+                                wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+                                constr, &signaller,
+                                lastbox,
+                                &totalNumberOfBondedInteractions, &bSumEkinhOld,
+                                (bGStat ? CGLO_GSTAT : 0)
+                                | (!EI_VV(ir->eI) || bRerunMD ? CGLO_ENERGY : 0)
+                                | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
+                                | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
+                                | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
+                                | CGLO_CONSTRAINT
+                                | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
+                                );
+                checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
+                                                top_global, top, state,
+                                                &shouldCheckNumberOfBondedInteractions);
+            }
         }
 
         /* #############  END CALC EKIN AND PRESSURE ################# */
@@ -1843,7 +1853,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         /* If it is time to reset counters, set a flag that remains
            true until counters actually get reset */
         if (step_rel == wcycle_get_reset_counters(wcycle) ||
-            gs.set[eglsRESETCOUNTERS] != 0)
+            signals[eglsRESETCOUNTERS].set != 0)
         {
             if (pme_loadbal_is_active(pme_loadbal))
             {
@@ -1872,9 +1882,9 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             /* Correct max_hours for the elapsed time */
             max_hours                -= elapsed_time/(60.0*60.0);
             /* If mdrun -maxh -resethway was active, it can only trigger once */
-            bResetCountersHalfMaxH    = FALSE; /* TODO move this to where gs.sig[eglsRESETCOUNTERS] is set */
+            bResetCountersHalfMaxH    = FALSE; /* TODO move this to where signals[eglsRESETCOUNTERS].sig is set */
             /* Reset can only happen once, so clear the triggering flag. */
-            gs.set[eglsRESETCOUNTERS] = 0;
+            signals[eglsRESETCOUNTERS].set = 0;
         }
 
         /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
@@ -1930,11 +1940,6 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         finish_swapcoords(ir->swap);
     }
 
-    if (membed != nullptr)
-    {
-        free_membed(membed);
-    }
-
     /* IMD cleanup, if bIMD is TRUE. */
     IMD_finalize(ir->bIMD, ir->imd);
 
diff --git a/patches/gromacs-2016-beta2.diff/src/programs/mdrun/md.cpp.preplumed b/patches/gromacs-2016.diff/src/programs/mdrun/md.cpp.preplumed
similarity index 92%
rename from patches/gromacs-2016-beta2.diff/src/programs/mdrun/md.cpp.preplumed
rename to patches/gromacs-2016.diff/src/programs/mdrun/md.cpp.preplumed
index 92fa911eb..2a49aba4a 100644
--- a/patches/gromacs-2016-beta2.diff/src/programs/mdrun/md.cpp.preplumed
+++ b/patches/gromacs-2016.diff/src/programs/mdrun/md.cpp.preplumed
@@ -76,13 +76,13 @@
 #include "gromacs/mdlib/mdebin.h"
 #include "gromacs/mdlib/mdoutf.h"
 #include "gromacs/mdlib/mdrun.h"
-#include "gromacs/mdlib/mdrun_signalling.h"
 #include "gromacs/mdlib/nb_verlet.h"
 #include "gromacs/mdlib/nbnxn_gpu_data_mgmt.h"
 #include "gromacs/mdlib/ns.h"
 #include "gromacs/mdlib/shellfc.h"
 #include "gromacs/mdlib/sighandler.h"
 #include "gromacs/mdlib/sim_util.h"
+#include "gromacs/mdlib/simulationsignal.h"
 #include "gromacs/mdlib/tgroup.h"
 #include "gromacs/mdlib/trajectory_writing.h"
 #include "gromacs/mdlib/update.h"
@@ -124,6 +124,8 @@
 #include "corewrap.h"
 #endif
 
+using gmx::SimulationSignaller;
+
 /*! \brief Check whether bonded interactions are missing, if appropriate
  *
  * \param[in]    fplog                                  Log file pointer
@@ -217,6 +219,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                   t_nrnb *nrnb, gmx_wallcycle_t wcycle,
                   gmx_edsam_t ed, t_forcerec *fr,
                   int repl_ex_nst, int repl_ex_nex, int repl_ex_seed,
+                  gmx_membed_t *membed,
                   real cpt_period, real max_hours,
                   int imdport,
                   unsigned long Flags,
@@ -229,7 +232,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     gmx_bool        bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
     gmx_bool        bNS, bNStList, bSimAnn, bStopCM, bRerunMD,
                     bFirstStep, startingFromCheckpoint, bInitStep, bLastStep = FALSE,
-                    bBornRadii;
+                    bBornRadii, bUsingEnsembleRestraints;
     gmx_bool          bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
     gmx_bool          do_ene, do_log, do_verbose, bRerunWarnNoV = TRUE,
                       bForceUpdate = FALSE, bCPT;
@@ -252,7 +255,6 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     gmx_global_stat_t gstat;
     gmx_update_t     *upd   = NULL;
     t_graph          *graph = NULL;
-    gmx_signalling_t  gs;
     gmx_groups_t     *groups;
     gmx_ekindata_t   *ekind;
     gmx_shellfc_t    *shellfc;
@@ -273,9 +275,8 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     int             **trotter_seq;
     char              sbuf[STEPSTRSIZE], sbuf2[STEPSTRSIZE];
     int               handled_stop_condition = gmx_stop_cond_none; /* compare to get_stop_condition*/
-    gmx_int64_t       multisim_nsteps        = -1;                 /* number of steps to do  before first multisim
-                                                                          simulation stops. If equal to zero, don't
-                                                                          communicate any more between multisims.*/
+
+
     /* PME load balancing data for GPU kernels */
     pme_load_balancing_t *pme_loadbal      = NULL;
     gmx_bool              bPMETune         = FALSE;
@@ -283,7 +284,6 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
 
     /* Interactive MD */
     gmx_bool          bIMDstep = FALSE;
-    gmx_membed_t     *membed   = NULL;
 
 #ifdef GMX_FAHCORE
     /* Temporary addition for FAHCORE checkpointing */
@@ -295,8 +295,13 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
        code. So we do that alongside the first global energy reduction
        after a new DD is made. These variables handle whether the
        check happens, and the result it returns. */
-    bool shouldCheckNumberOfBondedInteractions = false;
-    int  totalNumberOfBondedInteractions       = -1;
+    bool              shouldCheckNumberOfBondedInteractions = false;
+    int               totalNumberOfBondedInteractions       = -1;
+
+    SimulationSignals signals;
+    // Most global communnication stages don't propagate mdrun
+    // signals, and will use this object to achieve that.
+    SimulationSignaller nullSignaller(nullptr, nullptr, false, false);
 
     /* Check for special mdrun options */
     bRerunMD = (Flags & MD_RERUN);
@@ -335,18 +340,6 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     }
     groups = &top_global->groups;
 
-    if (opt2bSet("-membed", nfile, fnm))
-    {
-        if (MASTER(cr))
-        {
-            fprintf(stderr, "Initializing membed");
-        }
-        /* Note that membed cannot work in parallel because mtop is
-         * changed here. Fix this if we ever want to make it run with
-         * multiple ranks. */
-        membed = init_membed(fplog, nfile, fnm, top_global, ir, state_global, cr, &cpt_period);
-    }
-
     if (ir->eSwapCoords != eswapNO)
     {
         /* Initialize ion swapping code */
@@ -592,7 +585,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     bSumEkinhOld = FALSE;
     compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
                     NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
-                    constr, NULL, FALSE, state->box,
+                    constr, &nullSignaller, state->box,
                     &totalNumberOfBondedInteractions, &bSumEkinhOld, cglo_flags
                     | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0));
     checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
@@ -608,7 +601,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
 
         compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
                         NULL, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
-                        constr, NULL, FALSE, state->box,
+                        constr, &nullSignaller, state->box,
                         NULL, &bSumEkinhOld,
                         cglo_flags &~(CGLO_STOPCM | CGLO_PRESSURE));
     }
@@ -759,22 +752,46 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
     bSumEkinhOld     = FALSE;
     bExchanged       = FALSE;
     bNeedRepartition = FALSE;
+    // TODO This implementation of ensemble orientation restraints is nasty because
+    // a user can't just do multi-sim with single-sim orientation restraints.
+    bUsingEnsembleRestraints = (fcd->disres.nsystems > 1) || (cr->ms && fcd->orires.nr);
 
-    init_global_signals(&gs, cr, ir, repl_ex_nst);
+    {
+        // Replica exchange and ensemble restraints need all
+        // simulations to remain synchronized, so they need
+        // checkpoints and stop conditions to act on the same step, so
+        // the propagation of such signals must take place between
+        // simulations, not just within simulations.
+        bool checkpointIsLocal    = (repl_ex_nst <= 0) && !bUsingEnsembleRestraints;
+        bool stopConditionIsLocal = (repl_ex_nst <= 0) && !bUsingEnsembleRestraints;
+        bool resetCountersIsLocal = true;
+        signals[eglsCHKPT]         = SimulationSignal(checkpointIsLocal);
+        signals[eglsSTOPCOND]      = SimulationSignal(stopConditionIsLocal);
+        signals[eglsRESETCOUNTERS] = SimulationSignal(resetCountersIsLocal);
+    }
 
     step     = ir->init_step;
     step_rel = 0;
 
-    if (MULTISIM(cr) && (repl_ex_nst <= 0 ))
+    // TODO extract this to new multi-simulation module
+    if (MASTER(cr) && MULTISIM(cr) && (repl_ex_nst <= 0 ))
     {
-        /* check how many steps are left in other sims */
-        multisim_nsteps = get_multisim_nsteps(cr, ir->nsteps);
+        if (!multisim_int_all_are_equal(cr->ms, ir->nsteps))
+        {
+            md_print_info(cr, fplog,
+                          "Note: The number of steps is not consistent across multi simulations,\n"
+                          "but we are proceeding anyway!\n");
+        }
+        if (!multisim_int_all_are_equal(cr->ms, ir->init_step))
+        {
+            md_print_info(cr, fplog,
+                          "Note: The initial step is not consistent across multi simulations,\n"
+                          "but we are proceeding anyway!\n");
+        }
     }
 
-
     /* and stop now if we should */
-    bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps) ||
-                 ((multisim_nsteps >= 0) && (step_rel >= multisim_nsteps )));
+    bLastStep = (bLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
     while (!bLastStep)
     {
 
@@ -817,6 +834,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             t         = t0 + step*ir->delta_t;
         }
 
+        // TODO Refactor this, so that nstfep does not need a default value of zero
         if (ir->efep != efepNO || ir->bSimTemp)
         {
             /* find and set the current lambdas.  If rerunning, we either read in a state, or a lambda value,
@@ -902,7 +920,6 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         {
             /* for rerun MD always do Neighbour Searching */
             bNS      = (bFirstStep || ir->nstlist != 0);
-            bNStList = bNS;
         }
         else
         {
@@ -910,30 +927,9 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             bNS = (bFirstStep || bNStList || bExchanged || bNeedRepartition);
         }
 
-        /* check whether we should stop because another simulation has
-           stopped. */
-        if (MULTISIM(cr))
-        {
-            if ( (multisim_nsteps >= 0) &&  (step_rel >= multisim_nsteps)  &&
-                 (multisim_nsteps != ir->nsteps) )
-            {
-                if (bNS)
-                {
-                    if (MASTER(cr))
-                    {
-                        fprintf(stderr,
-                                "Stopping simulation %d because another one has finished\n",
-                                cr->ms->sim);
-                    }
-                    bLastStep         = TRUE;
-                    gs.sig[eglsCHKPT] = 1;
-                }
-            }
-        }
-
-        /* < 0 means stop after this step, > 0 means stop at next NS step */
-        if ( (gs.set[eglsSTOPCOND] < 0) ||
-             ( (gs.set[eglsSTOPCOND] > 0) && (bNStList || ir->nstlist == 0) ) )
+        /* < 0 means stop at next step, > 0 means stop at next NS step */
+        if ( (signals[eglsSTOPCOND].set < 0) ||
+             ( (signals[eglsSTOPCOND].set > 0 ) && ( bNS || ir->nstlist == 0)))
         {
             bLastStep = TRUE;
         }
@@ -1011,7 +1007,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             /* This may not be quite working correctly yet . . . . */
             compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
                             wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
-                            constr, NULL, FALSE, state->box,
+                            constr, &nullSignaller, state->box,
                             &totalNumberOfBondedInteractions, &bSumEkinhOld,
                             CGLO_GSTAT | CGLO_TEMPERATURE | CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS);
             checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
@@ -1025,12 +1021,12 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
          * or at the last step (but not when we do not want confout),
          * but never at the first step or with rerun.
          */
-        bCPT = (((gs.set[eglsCHKPT] && (bNS || ir->nstlist == 0)) ||
+        bCPT = (((signals[eglsCHKPT].set && (bNS || ir->nstlist == 0)) ||
                  (bLastStep && (Flags & MD_CONFOUT))) &&
                 step > ir->init_step && !bRerunMD);
         if (bCPT)
         {
-            gs.set[eglsCHKPT] = 0;
+            signals[eglsCHKPT].set = 0;
         }
 
         /* Determine the energy and pressure:
@@ -1168,7 +1164,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                 wallcycle_stop(wcycle, ewcUPDATE);
                 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
-                                constr, NULL, FALSE, state->box,
+                                constr, &nullSignaller, state->box,
                                 &totalNumberOfBondedInteractions, &bSumEkinhOld,
                                 (bGStat ? CGLO_GSTAT : 0)
                                 | CGLO_ENERGY
@@ -1216,7 +1212,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                     /* This may not be quite working correctly yet . . . . */
                     compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
                                     wcycle, enerd, NULL, NULL, NULL, NULL, mu_tot,
-                                    constr, NULL, FALSE, state->box,
+                                    constr, &nullSignaller, state->box,
                                     NULL, &bSumEkinhOld,
                                     CGLO_GSTAT | CGLO_TEMPERATURE);
                     wallcycle_start(wcycle, ewcUPDATE);
@@ -1295,18 +1291,18 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         {
             int nsteps_stop = -1;
 
-            /* this is just make gs.sig compatible with the hack
-               of sending signals around by MPI_Reduce with together with
+            /* this just makes signals[].sig compatible with the hack
+               of sending signals around by MPI_Reduce together with
                other floats */
             if (gmx_get_stop_condition() == gmx_stop_cond_next_ns)
             {
-                gs.sig[eglsSTOPCOND] = 1;
-                nsteps_stop          = std::max(ir->nstlist, 2*nstglobalcomm);
+                signals[eglsSTOPCOND].sig = 1;
+                nsteps_stop               = std::max(ir->nstlist, 2*nstglobalcomm);
             }
             if (gmx_get_stop_condition() == gmx_stop_cond_next)
             {
-                gs.sig[eglsSTOPCOND] = -1;
-                nsteps_stop          = nstglobalcomm + 1;
+                signals[eglsSTOPCOND].sig = -1;
+                nsteps_stop               = nstglobalcomm + 1;
             }
             if (fplog)
             {
@@ -1323,10 +1319,10 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         }
         else if (MASTER(cr) && (bNS || ir->nstlist <= 0) &&
                  (max_hours > 0 && elapsed_time > max_hours*60.0*60.0*0.99) &&
-                 gs.sig[eglsSTOPCOND] == 0 && gs.set[eglsSTOPCOND] == 0)
+                 signals[eglsSTOPCOND].sig == 0 && signals[eglsSTOPCOND].set == 0)
         {
             /* Signal to terminate the run */
-            gs.sig[eglsSTOPCOND] = 1;
+            signals[eglsSTOPCOND].sig = 1;
             if (fplog)
             {
                 fprintf(fplog, "\nStep %s: Run time exceeded %.3f hours, will terminate the run\n", gmx_step_str(step, sbuf), max_hours*0.99);
@@ -1338,7 +1334,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             elapsed_time > max_hours*60.0*60.0*0.495)
         {
             /* Set flag that will communicate the signal to all ranks in the simulation */
-            gs.sig[eglsRESETCOUNTERS] = 1;
+            signals[eglsRESETCOUNTERS].sig = 1;
         }
 
         /* In parallel we only have to check for checkpointing in steps
@@ -1349,9 +1345,9 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                            cpt_period >= 0 &&
                            (cpt_period == 0 ||
                             elapsed_time >= nchkpt*cpt_period*60.0)) &&
-            gs.set[eglsCHKPT] == 0)
+            signals[eglsCHKPT].set == 0)
         {
-            gs.sig[eglsCHKPT] = 1;
+            signals[eglsCHKPT].sig = 1;
         }
 
         /* #########   START SECOND UPDATE STEP ################# */
@@ -1441,7 +1437,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
                 /* just compute the kinetic energy at the half step to perform a trotter step */
                 compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
                                 wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
-                                constr, NULL, FALSE, lastbox,
+                                constr, &nullSignaller, lastbox,
                                 NULL, &bSumEkinhOld,
                                 (bGStat ? CGLO_GSTAT : 0) | CGLO_TEMPERATURE
                                 );
@@ -1517,26 +1513,40 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
          * non-communication steps, but we need to calculate
          * the kinetic energy one step before communication.
          */
-        if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)))
         {
-            compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
-                            wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
-                            constr, &gs,
-                            (step_rel % gs.nstms == 0) &&
-                            (multisim_nsteps < 0 || (step_rel < multisim_nsteps)),
-                            lastbox,
-                            &totalNumberOfBondedInteractions, &bSumEkinhOld,
-                            (bGStat ? CGLO_GSTAT : 0)
-                            | (!EI_VV(ir->eI) || bRerunMD ? CGLO_ENERGY : 0)
-                            | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
-                            | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
-                            | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
-                            | CGLO_CONSTRAINT
-                            | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
-                            );
-            checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
-                                            top_global, top, state,
-                                            &shouldCheckNumberOfBondedInteractions);
+            // Organize to do inter-simulation signalling on steps if
+            // and when algorithms require it.
+            bool doInterSimSignal = (!bFirstStep && bDoReplEx) || bUsingEnsembleRestraints;
+
+            if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)) || doInterSimSignal)
+            {
+                // Since we're already communicating at this step, we
+                // can propagate intra-simulation signals. Note that
+                // check_nstglobalcomm has the responsibility for
+                // choosing the value of nstglobalcomm that is one way
+                // bGStat becomes true, so we can't get into a
+                // situation where e.g. checkpointing can't be
+                // signalled.
+                bool                doIntraSimSignal = true;
+                SimulationSignaller signaller(&signals, cr, doInterSimSignal, doIntraSimSignal);
+
+                compute_globals(fplog, gstat, cr, ir, fr, ekind, state, mdatoms, nrnb, vcm,
+                                wcycle, enerd, force_vir, shake_vir, total_vir, pres, mu_tot,
+                                constr, &signaller,
+                                lastbox,
+                                &totalNumberOfBondedInteractions, &bSumEkinhOld,
+                                (bGStat ? CGLO_GSTAT : 0)
+                                | (!EI_VV(ir->eI) || bRerunMD ? CGLO_ENERGY : 0)
+                                | (!EI_VV(ir->eI) && bStopCM ? CGLO_STOPCM : 0)
+                                | (!EI_VV(ir->eI) ? CGLO_TEMPERATURE : 0)
+                                | (!EI_VV(ir->eI) || bRerunMD ? CGLO_PRESSURE : 0)
+                                | CGLO_CONSTRAINT
+                                | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS : 0)
+                                );
+                checkNumberOfBondedInteractions(fplog, cr, totalNumberOfBondedInteractions,
+                                                top_global, top, state,
+                                                &shouldCheckNumberOfBondedInteractions);
+            }
         }
 
         /* #############  END CALC EKIN AND PRESSURE ################# */
@@ -1740,7 +1750,7 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         /* If it is time to reset counters, set a flag that remains
            true until counters actually get reset */
         if (step_rel == wcycle_get_reset_counters(wcycle) ||
-            gs.set[eglsRESETCOUNTERS] != 0)
+            signals[eglsRESETCOUNTERS].set != 0)
         {
             if (pme_loadbal_is_active(pme_loadbal))
             {
@@ -1769,9 +1779,9 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
             /* Correct max_hours for the elapsed time */
             max_hours                -= elapsed_time/(60.0*60.0);
             /* If mdrun -maxh -resethway was active, it can only trigger once */
-            bResetCountersHalfMaxH    = FALSE; /* TODO move this to where gs.sig[eglsRESETCOUNTERS] is set */
+            bResetCountersHalfMaxH    = FALSE; /* TODO move this to where signals[eglsRESETCOUNTERS].sig is set */
             /* Reset can only happen once, so clear the triggering flag. */
-            gs.set[eglsRESETCOUNTERS] = 0;
+            signals[eglsRESETCOUNTERS].set = 0;
         }
 
         /* If bIMD is TRUE, the master updates the IMD energy record and sends positions to VMD client */
@@ -1827,11 +1837,6 @@ double gmx::do_md(FILE *fplog, t_commrec *cr, int nfile, const t_filenm fnm[],
         finish_swapcoords(ir->swap);
     }
 
-    if (membed != nullptr)
-    {
-        free_membed(membed);
-    }
-
     /* IMD cleanup, if bIMD is TRUE. */
     IMD_finalize(ir->bIMD, ir->imd);
 
diff --git a/patches/gromacs-2016-beta2.diff/src/programs/mdrun/mdrun.cpp b/patches/gromacs-2016.diff/src/programs/mdrun/mdrun.cpp
similarity index 97%
rename from patches/gromacs-2016-beta2.diff/src/programs/mdrun/mdrun.cpp
rename to patches/gromacs-2016.diff/src/programs/mdrun/mdrun.cpp
index 3efc27c46..2c3906175 100644
--- a/patches/gromacs-2016-beta2.diff/src/programs/mdrun/mdrun.cpp
+++ b/patches/gromacs-2016.diff/src/programs/mdrun/mdrun.cpp
@@ -489,16 +489,30 @@ int gmx_mdrun(int argc, char *argv[])
     if (nmultisim >= 1)
     {
 #if !GMX_THREAD_MPI
-        gmx_bool bParFn = (multidir == NULL);
-        init_multisystem(cr, nmultisim, multidir, NFILE, fnm, bParFn);
+        init_multisystem(cr, nmultisim, multidir, NFILE, fnm);
 #else
         gmx_fatal(FARGS, "mdrun -multi or -multidir are not supported with the thread-MPI library. "
                   "Please compile GROMACS with a proper external MPI library.");
 #endif
     }
 
-    handleRestart(cr, bTryToAppendFiles, NFILE, fnm,
-                  &bDoAppendFiles, &bStartFromCpt);
+    if (!opt2bSet("-cpi", NFILE, fnm))
+    {
+        // If we are not starting from a checkpoint we never allow files to be appended
+        // to, since that has caused a ton of strange behaviour and bugs in the past.
+        if (opt2parg_bSet("-append", asize(pa), pa))
+        {
+            // If the user explicitly used the -append option, explain that it is not possible.
+            gmx_fatal(FARGS, "GROMACS can only append to files when restarting from a checkpoint.");
+        }
+        else
+        {
+            // If the user did not say anything explicit, just disable appending.
+            bTryToAppendFiles = FALSE;
+        }
+    }
+
+    handleRestart(cr, bTryToAppendFiles, NFILE, fnm, &bDoAppendFiles, &bStartFromCpt);
 
     Flags = opt2bSet("-rerun", NFILE, fnm) ? MD_RERUN : 0;
     Flags = Flags | (bDDBondCheck  ? MD_DDBONDCHECK  : 0);
@@ -507,7 +521,7 @@ int gmx_mdrun(int argc, char *argv[])
     Flags = Flags | (bConfout      ? MD_CONFOUT      : 0);
     Flags = Flags | (bRerunVSite   ? MD_RERUN_VSITE  : 0);
     Flags = Flags | (bReproducible ? MD_REPRODUCIBLE : 0);
-    Flags = Flags | (bDoAppendFiles ? MD_APPENDFILES  : 0);
+    Flags = Flags | (bDoAppendFiles  ? MD_APPENDFILES  : 0);
     Flags = Flags | (opt2parg_bSet("-append", asize(pa), pa) ? MD_APPENDFILESSET : 0);
     Flags = Flags | (bKeepAndNumCPT ? MD_KEEPANDNUMCPT : 0);
     Flags = Flags | (bStartFromCpt ? MD_STARTFROMCPT : 0);
diff --git a/patches/gromacs-2016-beta2.diff/src/programs/mdrun/mdrun.cpp.preplumed b/patches/gromacs-2016.diff/src/programs/mdrun/mdrun.cpp.preplumed
similarity index 96%
rename from patches/gromacs-2016-beta2.diff/src/programs/mdrun/mdrun.cpp.preplumed
rename to patches/gromacs-2016.diff/src/programs/mdrun/mdrun.cpp.preplumed
index e56c4a145..7829a823a 100644
--- a/patches/gromacs-2016-beta2.diff/src/programs/mdrun/mdrun.cpp.preplumed
+++ b/patches/gromacs-2016.diff/src/programs/mdrun/mdrun.cpp.preplumed
@@ -481,16 +481,30 @@ int gmx_mdrun(int argc, char *argv[])
     if (nmultisim >= 1)
     {
 #if !GMX_THREAD_MPI
-        gmx_bool bParFn = (multidir == NULL);
-        init_multisystem(cr, nmultisim, multidir, NFILE, fnm, bParFn);
+        init_multisystem(cr, nmultisim, multidir, NFILE, fnm);
 #else
         gmx_fatal(FARGS, "mdrun -multi or -multidir are not supported with the thread-MPI library. "
                   "Please compile GROMACS with a proper external MPI library.");
 #endif
     }
 
-    handleRestart(cr, bTryToAppendFiles, NFILE, fnm,
-                  &bDoAppendFiles, &bStartFromCpt);
+    if (!opt2bSet("-cpi", NFILE, fnm))
+    {
+        // If we are not starting from a checkpoint we never allow files to be appended
+        // to, since that has caused a ton of strange behaviour and bugs in the past.
+        if (opt2parg_bSet("-append", asize(pa), pa))
+        {
+            // If the user explicitly used the -append option, explain that it is not possible.
+            gmx_fatal(FARGS, "GROMACS can only append to files when restarting from a checkpoint.");
+        }
+        else
+        {
+            // If the user did not say anything explicit, just disable appending.
+            bTryToAppendFiles = FALSE;
+        }
+    }
+
+    handleRestart(cr, bTryToAppendFiles, NFILE, fnm, &bDoAppendFiles, &bStartFromCpt);
 
     Flags = opt2bSet("-rerun", NFILE, fnm) ? MD_RERUN : 0;
     Flags = Flags | (bDDBondCheck  ? MD_DDBONDCHECK  : 0);
@@ -499,7 +513,7 @@ int gmx_mdrun(int argc, char *argv[])
     Flags = Flags | (bConfout      ? MD_CONFOUT      : 0);
     Flags = Flags | (bRerunVSite   ? MD_RERUN_VSITE  : 0);
     Flags = Flags | (bReproducible ? MD_REPRODUCIBLE : 0);
-    Flags = Flags | (bDoAppendFiles ? MD_APPENDFILES  : 0);
+    Flags = Flags | (bDoAppendFiles  ? MD_APPENDFILES  : 0);
     Flags = Flags | (opt2parg_bSet("-append", asize(pa), pa) ? MD_APPENDFILESSET : 0);
     Flags = Flags | (bKeepAndNumCPT ? MD_KEEPANDNUMCPT : 0);
     Flags = Flags | (bStartFromCpt ? MD_STARTFROMCPT : 0);
diff --git a/patches/gromacs-2016-beta2.diff/src/programs/mdrun/repl_ex.cpp b/patches/gromacs-2016.diff/src/programs/mdrun/repl_ex.cpp
similarity index 100%
rename from patches/gromacs-2016-beta2.diff/src/programs/mdrun/repl_ex.cpp
rename to patches/gromacs-2016.diff/src/programs/mdrun/repl_ex.cpp
diff --git a/patches/gromacs-2016-beta2.diff/src/programs/mdrun/repl_ex.cpp.preplumed b/patches/gromacs-2016.diff/src/programs/mdrun/repl_ex.cpp.preplumed
similarity index 100%
rename from patches/gromacs-2016-beta2.diff/src/programs/mdrun/repl_ex.cpp.preplumed
rename to patches/gromacs-2016.diff/src/programs/mdrun/repl_ex.cpp.preplumed
diff --git a/patches/gromacs-2016-beta2.diff/src/programs/mdrun/runner.cpp b/patches/gromacs-2016.diff/src/programs/mdrun/runner.cpp
similarity index 98%
rename from patches/gromacs-2016-beta2.diff/src/programs/mdrun/runner.cpp
rename to patches/gromacs-2016.diff/src/programs/mdrun/runner.cpp
index b8550c39e..f3e0e4ece 100644
--- a/patches/gromacs-2016-beta2.diff/src/programs/mdrun/runner.cpp
+++ b/patches/gromacs-2016.diff/src/programs/mdrun/runner.cpp
@@ -110,6 +110,7 @@
 
 #include "deform.h"
 #include "md.h"
+#include "membed.h"
 #include "repl_ex.h"
 #include "resource-division.h"
 
@@ -715,6 +716,7 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
     gmx_int64_t               reset_counters;
     gmx_edsam_t               ed           = NULL;
     int                       nthreads_pme = 1;
+    gmx_membed_t *            membed       = NULL;
     gmx_hw_info_t            *hwinfo       = NULL;
     /* The master rank decides early on bUseGPU and broadcasts this later */
     gmx_bool                  bUseGPU            = FALSE;
@@ -729,6 +731,7 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
         fplog = NULL;
     }
 
+    bool doMembed = opt2bSet("-membed", nfile, fnm);
     bRerunMD     = (Flags & MD_RERUN);
     bForceUseGPU = (strncmp(nbpu_opt, "gpu", 3) == 0);
     bTryUseGPU   = (strncmp(nbpu_opt, "auto", 4) == 0) || bForceUseGPU;
@@ -836,7 +839,8 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
         hw_opt->nthreads_tmpi = get_nthreads_mpi(hwinfo,
                                                  hw_opt,
                                                  inputrec, mtop,
-                                                 cr, fplog, bUseGPU);
+                                                 cr, fplog, bUseGPU,
+                                                 doMembed);
 
         if (hw_opt->nthreads_tmpi > 1)
         {
@@ -1148,6 +1152,19 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
         wcycle_set_reset_counters(wcycle, reset_counters);
     }
 
+    // Membrane embedding must be initialized before we call init_forcerec()
+    if (doMembed)
+    {
+        if (MASTER(cr))
+        {
+            fprintf(stderr, "Initializing membed");
+        }
+        /* Note that membed cannot work in parallel because mtop is
+         * changed here. Fix this if we ever want to make it run with
+         * multiple ranks. */
+        membed = init_membed(fplog, nfile, fnm, mtop, inputrec, state, cr, &cpt_period);
+    }
+
     snew(nrnb, 1);
     if (cr->duty & DUTY_PP)
     {
@@ -1323,6 +1340,7 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
                                      fcd, state,
                                      mdatoms, nrnb, wcycle, ed, fr,
                                      repl_ex_nst, repl_ex_nex, repl_ex_seed,
+                                     membed,
                                      cpt_period, max_hours,
                                      imdport,
                                      Flags,
@@ -1361,6 +1379,11 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
     /* Free GPU memory and context */
     free_gpu_resources(fr, cr, &hwinfo->gpu_info, fr ? fr->gpu_opt : NULL);
 
+    if (doMembed)
+    {
+        free_membed(membed);
+    }
+
     gmx_hardware_info_free(hwinfo);
 
     /* Does what it says */
diff --git a/patches/gromacs-2016-beta2.diff/src/programs/mdrun/runner.cpp.preplumed b/patches/gromacs-2016.diff/src/programs/mdrun/runner.cpp.preplumed
similarity index 98%
rename from patches/gromacs-2016-beta2.diff/src/programs/mdrun/runner.cpp.preplumed
rename to patches/gromacs-2016.diff/src/programs/mdrun/runner.cpp.preplumed
index ddc5895c9..560b923ba 100644
--- a/patches/gromacs-2016-beta2.diff/src/programs/mdrun/runner.cpp.preplumed
+++ b/patches/gromacs-2016.diff/src/programs/mdrun/runner.cpp.preplumed
@@ -110,6 +110,7 @@
 
 #include "deform.h"
 #include "md.h"
+#include "membed.h"
 #include "repl_ex.h"
 #include "resource-division.h"
 
@@ -709,6 +710,7 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
     gmx_int64_t               reset_counters;
     gmx_edsam_t               ed           = NULL;
     int                       nthreads_pme = 1;
+    gmx_membed_t *            membed       = NULL;
     gmx_hw_info_t            *hwinfo       = NULL;
     /* The master rank decides early on bUseGPU and broadcasts this later */
     gmx_bool                  bUseGPU            = FALSE;
@@ -723,6 +725,7 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
         fplog = NULL;
     }
 
+    bool doMembed = opt2bSet("-membed", nfile, fnm);
     bRerunMD     = (Flags & MD_RERUN);
     bForceUseGPU = (strncmp(nbpu_opt, "gpu", 3) == 0);
     bTryUseGPU   = (strncmp(nbpu_opt, "auto", 4) == 0) || bForceUseGPU;
@@ -830,7 +833,8 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
         hw_opt->nthreads_tmpi = get_nthreads_mpi(hwinfo,
                                                  hw_opt,
                                                  inputrec, mtop,
-                                                 cr, fplog, bUseGPU);
+                                                 cr, fplog, bUseGPU,
+                                                 doMembed);
 
         if (hw_opt->nthreads_tmpi > 1)
         {
@@ -1142,6 +1146,19 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
         wcycle_set_reset_counters(wcycle, reset_counters);
     }
 
+    // Membrane embedding must be initialized before we call init_forcerec()
+    if (doMembed)
+    {
+        if (MASTER(cr))
+        {
+            fprintf(stderr, "Initializing membed");
+        }
+        /* Note that membed cannot work in parallel because mtop is
+         * changed here. Fix this if we ever want to make it run with
+         * multiple ranks. */
+        membed = init_membed(fplog, nfile, fnm, mtop, inputrec, state, cr, &cpt_period);
+    }
+
     snew(nrnb, 1);
     if (cr->duty & DUTY_PP)
     {
@@ -1317,6 +1334,7 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
                                      fcd, state,
                                      mdatoms, nrnb, wcycle, ed, fr,
                                      repl_ex_nst, repl_ex_nex, repl_ex_seed,
+                                     membed,
                                      cpt_period, max_hours,
                                      imdport,
                                      Flags,
@@ -1355,6 +1373,11 @@ int mdrunner(gmx_hw_opt_t *hw_opt,
     /* Free GPU memory and context */
     free_gpu_resources(fr, cr, &hwinfo->gpu_info, fr ? fr->gpu_opt : NULL);
 
+    if (doMembed)
+    {
+        free_membed(membed);
+    }
+
     gmx_hardware_info_free(hwinfo);
 
     /* Does what it says */
-- 
GitLab