diff --git a/cmake/gmxManageMimic.cmake b/cmake/gmxManageMimic.cmake
index c0c98980bf..b5d28965b3 100644
--- a/cmake/gmxManageMimic.cmake
+++ b/cmake/gmxManageMimic.cmake
@@ -33,7 +33,7 @@
 # the research papers on the package. Check out http://www.gromacs.org.
 
 if(GMX_MIMIC)
-    find_package(MimicCommLib REQUIRED)
-    include_directories(SYSTEM ${MimicCommLib_INCLUDE_DIRS})
-    list(APPEND GMX_COMMON_LIBRARIES Upstream::mimiccomm)
+    find_package(MiMiCCommLib REQUIRED CXX)
+    include_directories(SYSTEM ${MiMiCCommLib_INCLUDE_DIRS})
+    list(APPEND GMX_COMMON_LIBRARIES ${MiMiCCommLib_LIBRARIES})
 endif()
diff --git a/src/gromacs/gmxlib/network.cpp b/src/gromacs/gmxlib/network.cpp
index 66023edac8..6b83ce7c03 100644
--- a/src/gromacs/gmxlib/network.cpp
+++ b/src/gromacs/gmxlib/network.cpp
@@ -114,12 +114,12 @@ void done_commrec(t_commrec* cr)
     // structure of the commrec and domdec initialization code makes
     // it hard to avoid both leaks and double frees.
     bool mySimIsMyGroup = (cr->mpi_comm_mysim == cr->mpi_comm_mygroup);
-    if (cr->mpi_comm_mysim != MPI_COMM_NULL && cr->mpi_comm_mysim != MPI_COMM_WORLD)
+    if (cr->mpi_comm_mysim != MPI_COMM_NULL && cr->mpi_comm_mysim != MPI_WORLD_STUB)
     {
         // TODO see above
         // MPI_Comm_free(&cr->mpi_comm_mysim);
     }
-    if (!mySimIsMyGroup && cr->mpi_comm_mygroup != MPI_COMM_NULL && cr->mpi_comm_mygroup != MPI_COMM_WORLD)
+    if (!mySimIsMyGroup && cr->mpi_comm_mygroup != MPI_COMM_NULL && cr->mpi_comm_mygroup != MPI_WORLD_STUB)
     {
         // TODO see above
         // MPI_Comm_free(&cr->mpi_comm_mygroup);
@@ -485,7 +485,7 @@ void gmx_fatal_collective(int                    f_errno,
 #if GMX_MPI
     int result;
     /* Check if we are calling on all processes in MPI_COMM_WORLD */
-    MPI_Comm_compare(comm, MPI_COMM_WORLD, &result);
+    MPI_Comm_compare(comm, MPI_WORLD_STUB, &result);
     /* Any result except MPI_UNEQUAL allows us to call MPI_Finalize */
     bFinalize = (result != MPI_UNEQUAL);
 #else
diff --git a/src/gromacs/hardware/detecthardware.cpp b/src/gromacs/hardware/detecthardware.cpp
index d808249678..8b2aa4cd3b 100644
--- a/src/gromacs/hardware/detecthardware.cpp
+++ b/src/gromacs/hardware/detecthardware.cpp
@@ -252,7 +252,7 @@ static void gmx_collect_hardware_mpi(const gmx::CpuInfo&             cpuInfo,
         }
 
         MPI_Allreduce(countsLocal.data(), countsReduced.data(), countsLocal.size(), MPI_INT,
-                      MPI_SUM, MPI_COMM_WORLD);
+                      MPI_SUM, physicalNodeComm.comm_);
     }
 
     constexpr int                   numElementsMax = 11;
@@ -275,7 +275,7 @@ static void gmx_collect_hardware_mpi(const gmx::CpuInfo&             cpuInfo,
         maxMinLocal[10] = (cpuIsAmdZen1 ? 1 : 0);
 
         MPI_Allreduce(maxMinLocal.data(), maxMinReduced.data(), maxMinLocal.size(), MPI_INT,
-                      MPI_MAX, MPI_COMM_WORLD);
+                      MPI_MAX, physicalNodeComm.comm_);
     }
 
     hardwareInfo->nphysicalnode       = countsReduced[0];
diff --git a/src/gromacs/hardware/printhardware.cpp b/src/gromacs/hardware/printhardware.cpp
index 0720a12bbe..398a03b8a2 100644
--- a/src/gromacs/hardware/printhardware.cpp
+++ b/src/gromacs/hardware/printhardware.cpp
@@ -213,7 +213,7 @@ static std::string detected_hardware_string(const gmx_hw_info_t* hwinfo, bool bF
 
     gmx_gethostname(host, STRLEN);
 
-    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+    MPI_Comm_rank(MPI_WORLD_STUB, &rank);
 
     // TODO Use a wrapper around MPI_Get_processor_name instead.
     s += gmx::formatString("Hardware detected on host %s (the node of MPI rank %d):\n", host, rank);
diff --git a/src/gromacs/mdrun/mimic.cpp b/src/gromacs/mdrun/mimic.cpp
index e37d921c5a..e49649281a 100644
--- a/src/gromacs/mdrun/mimic.cpp
+++ b/src/gromacs/mdrun/mimic.cpp
@@ -50,6 +50,8 @@
 #include <algorithm>
 #include <memory>
 
+#include "Requests.h"
+
 #include "gromacs/applied_forces/awh/awh.h"
 #include "gromacs/commandline/filenm.h"
 #include "gromacs/domdec/collect.h"
@@ -209,12 +211,11 @@ void gmx::LegacySimulator::do_mimic()
     int        nstglobalcomm = 1;
     const bool bNS           = true;
 
+    MimicCommunicator mimicCommunicator{};
     if (MASTER(cr))
     {
-        MimicCommunicator::init();
-        auto nonConstGlobalTopology = const_cast<gmx_mtop_t*>(top_global);
-        MimicCommunicator::sendInitData(nonConstGlobalTopology, state_global->x);
-        ir->nsteps = MimicCommunicator::getStepNumber();
+        mimicCommunicator.init();
+        ir->nsteps=2;
     }
 
     ir->nstxout_compressed         = 0;
@@ -266,7 +267,10 @@ void gmx::LegacySimulator::do_mimic()
                             imdSession, pull_work, state, &f, mdAtoms, &top, fr, vsite, constr,
                             nrnb, nullptr, FALSE);
         shouldCheckNumberOfBondedInteractions = true;
-        gmx_bcast(sizeof(ir->nsteps), &ir->nsteps, cr->mpi_comm_mygroup);
+        if (PAR(cr))
+        {
+            gmx_bcast(sizeof(ir->nsteps), &ir->nsteps, cr->mpi_comm_mygroup);
+        }
     }
     else
     {
@@ -361,122 +365,124 @@ void gmx::LegacySimulator::do_mimic()
     isLastStep = (isLastStep || (ir->nsteps >= 0 && step_rel > ir->nsteps));
     while (!isLastStep)
     {
-        isLastStep = (isLastStep || (ir->nsteps >= 0 && step_rel == ir->nsteps));
-        wallcycle_start(wcycle, ewcSTEP);
-
-        t = step;
-
+        int command = -1;
         if (MASTER(cr))
         {
-            MimicCommunicator::getCoords(&state_global->x, state_global->natoms);
+            command = mimicCommunicator.getCommand();
         }
 
-        if (ir->efep != efepNO)
+        if (PAR(cr))
         {
-            state->lambda = currentLambdas(step, *(ir->fepvals), state_global->fep_state);
+            gmx_barrier(cr->mpi_comm_mygroup);
+            gmx_bcast(sizeof(int), &command, cr->mpi_comm_mygroup);
         }
 
-        if (MASTER(cr))
+        t = step;
+
+        auto* nonConstGlobalTopology = const_cast<gmx_mtop_t*>(top_global);
+        if (command == MCL_EXIT)
         {
-            const bool constructVsites = ((vsite != nullptr) && mdrunOptions.rerunConstructVsites);
-            if (constructVsites && DOMAINDECOMP(cr))
+            isLastStep = true;
+        }
+        else if (command == MCL_SEND_ATOM_SPECIES)
+        {
+            if (MASTER(cr))
             {
-                gmx_fatal(FARGS,
-                          "Vsite recalculation with -rerun is not implemented with domain "
-                          "decomposition, "
-                          "use a single rank");
+                mimicCommunicator.sendAtomSpecies(nonConstGlobalTopology);
             }
-            if (constructVsites)
+        }
+        else if (command == MCL_SEND_NUM_FRAGMENTS)
+        {
+            if (MASTER(cr))
             {
-                wallcycle_start(wcycle, ewcVSITECONSTR);
-                vsite->construct(state->x, ir->delta_t, state->v, state->box);
-                wallcycle_stop(wcycle, ewcVSITECONSTR);
+                mimicCommunicator.sendFragmentCount(nonConstGlobalTopology);
             }
         }
-
-        if (DOMAINDECOMP(cr))
+        else if (command == MCL_SEND_BOND_ATOMS)
         {
-            /* Repartition the domain decomposition */
-            const bool bMasterState = true;
-            dd_partition_system(fplog, mdlog, step, cr, bMasterState, nstglobalcomm, state_global,
-                                *top_global, ir, imdSession, pull_work, state, &f, mdAtoms, &top,
-                                fr, vsite, constr, nrnb, wcycle, mdrunOptions.verbose);
-            shouldCheckNumberOfBondedInteractions = true;
+            if (MASTER(cr))
+            {
+                mimicCommunicator.sendBondAtoms(nonConstGlobalTopology);
+            }
         }
-
-        if (MASTER(cr))
+        else if (command == MCL_SEND_BOND_LENGTHS)
         {
-            EnergyOutput::printHeader(fplog, step, t); /* can we improve the information printed here? */
+            if (MASTER(cr))
+            {
+                mimicCommunicator.sendBondLengths(nonConstGlobalTopology);
+            }
         }
-
-        if (ir->efep != efepNO)
+        else if (command == MCL_SEND_ATOM_ELEMENTS)
         {
-            update_mdatoms(mdatoms, state->lambda[efptMASS]);
+            if (MASTER(cr))
+            {
+                mimicCommunicator.sendAtomElements(nonConstGlobalTopology);
+            }
         }
-
-        force_flags = (GMX_FORCE_STATECHANGED | GMX_FORCE_DYNAMICBOX | GMX_FORCE_ALLFORCES
-                       | GMX_FORCE_VIRIAL | // TODO: Get rid of this once #2649 is solved
-                       GMX_FORCE_ENERGY | (doFreeEnergyPerturbation ? GMX_FORCE_DHDL : 0));
-
-        if (shellfc)
+        else if (command == MCL_SEND_NUM_ATOMS_IN_FRAGMENTS)
         {
-            /* Now is the time to relax the shells */
-            relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, enforcedRotation, step, ir,
-                                imdSession, pull_work, bNS, force_flags, &top, constr, enerd,
-                                state->natoms, state->x.arrayRefWithPadding(),
-                                state->v.arrayRefWithPadding(), state->box, state->lambda,
-                                &state->hist, &f.view(), force_vir, mdatoms, nrnb, wcycle, shellfc,
-                                fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler);
+            if (MASTER(cr))
+            {
+                mimicCommunicator.sendFragmentAtomNum(nonConstGlobalTopology);
+            }
         }
-        else
+        else if (command == MCL_SEND_ATOMS_IN_FRAGMENTS)
         {
-            /* The coordinates (x) are shifted (to get whole molecules)
-             * in do_force.
-             * This is parallellized as well, and does communication too.
-             * Check comments in sim_util.c
-             */
-            Awh*       awh = nullptr;
-            gmx_edsam* ed  = nullptr;
-            do_force(fplog, cr, ms, ir, awh, enforcedRotation, imdSession, pull_work, step, nrnb,
-                     wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist,
-                     &f.view(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork,
-                     vsite, mu_tot, t, ed, GMX_FORCE_NS | force_flags, ddBalanceRegionHandler);
+            if (MASTER(cr))
+            {
+                mimicCommunicator.sendFragmentAtomIds(nonConstGlobalTopology);
+            }
         }
-
-        /* Now we have the energies and forces corresponding to the
-         * coordinates at time t.
-         */
+        else if (command == MCL_SEND_ATOM_MASSES)
         {
-            const bool isCheckpointingStep = false;
-            const bool doRerun             = false;
-            const bool bSumEkinhOld        = false;
-            do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t, ir, state,
-                                     state_global, observablesHistory, top_global, fr, outf,
-                                     energyOutput, ekind, f.view().force(), isCheckpointingStep,
-                                     doRerun, isLastStep, mdrunOptions.writeConfout, bSumEkinhOld);
+            if (MASTER(cr))
+            {
+                mimicCommunicator.sendAtomMasses(nonConstGlobalTopology);
+            }
         }
-
-        stopHandler->setSignal();
-
+        else if (command == MCL_SEND_ATOM_MULTIPOLES)
+        {
+            if (MASTER(cr))
+            {
+                mimicCommunicator.sendCharges(nonConstGlobalTopology);
+            }
+        }
+        else if (command == MCL_SEND_ATOM_COORDS)
         {
-            const bool          doInterSimSignal = false;
-            const bool          doIntraSimSignal = true;
-            bool                bSumEkinhOld     = false;
-            t_vcm*              vcm              = nullptr;
-            SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
-
-            compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
-                            makeConstArrayRef(state->v), state->box, mdatoms, nrnb, vcm, wcycle,
-                            enerd, nullptr, nullptr, nullptr, nullptr, constr, &signaller,
-                            state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
-                            CGLO_GSTAT | CGLO_ENERGY
-                                    | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
-                                                                             : 0));
-            checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions, top_global,
-                                            &top, makeConstArrayRef(state->x), state->box,
-                                            &shouldCheckNumberOfBondedInteractions);
+            if (MASTER(cr))
+            {
+                mimicCommunicator.sendCoords(state_global->x);
+            }
         }
-
+        else if (command == MCL_SEND_NUM_ATOMS)
+        {
+            if (MASTER(cr))
+            {
+                mimicCommunicator.sendAtomNum(nonConstGlobalTopology);
+            }
+        }
+        else if (command == MCL_SEND_NUM_ATOM_TYPES)
+        {
+            if (MASTER(cr))
+            {
+                mimicCommunicator.sendSpeciesNum(nonConstGlobalTopology);
+            }
+        }
+        else if (command == MCL_SEND_NUM_BONDS)
+        {
+            if (MASTER(cr))
+            {
+                mimicCommunicator.sendBondsNumber(nonConstGlobalTopology);
+            }
+        }
+        else if (command == MCL_SEND_NUM_ANGLES)
+        {
+            if (MASTER(cr))
+            {
+                mimicCommunicator.sendAnglesNumber(nonConstGlobalTopology);
+            }
+        }
+        else if (command == MCL_SEND_ATOM_FORCES)
         {
             gmx::HostVector<gmx::RVec>     fglobal(top_global->natoms);
             gmx::ArrayRef<gmx::RVec>       ftemp;
@@ -494,73 +500,206 @@ void gmx::LegacySimulator::do_mimic()
 
             if (MASTER(cr))
             {
-                MimicCommunicator::sendEnergies(enerd->term[F_EPOT]);
-                MimicCommunicator::sendForces(ftemp, state_global->natoms);
+                mimicCommunicator.sendForces(ftemp, state_global->natoms);
             }
         }
-
-
-        /* Note: this is OK, but there are some numerical precision issues with using the convergence of
-           the virial that should probably be addressed eventually. state->veta has better properies,
-           but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
-           generate the new shake_vir, but test the veta value for convergence.  This will take some thought. */
-
-        if (ir->efep != efepNO)
+        else if (command == MCL_SEND_ENERGY)
         {
-            /* Sum up the foreign energy and dhdl terms for md and sd.
-               Currently done every step so that dhdl is correct in the .edr */
-            accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
+            if (MASTER(cr))
+            {
+                mimicCommunicator.sendEnergies(enerd->term[F_EPOT]);
+            }
         }
-
-        /* Output stuff */
-        if (MASTER(cr))
+        else if (command == MCL_SEND_MULTIPOLE_ORDER)
         {
-            const bool bCalcEnerStep = true;
-            energyOutput.addDataAtEnergyStep(
-                    doFreeEnergyPerturbation, bCalcEnerStep, t, mdatoms->tmass, enerd, ir->fepvals,
-                    ir->expandedvals, state->box,
-                    PTCouplingArrays({ state->boxv, state->nosehoover_xi, state->nosehoover_vxi,
-                                       state->nhpres_xi, state->nhpres_vxi }),
-                    state->fep_state, shake_vir, force_vir, total_vir, pres, ekind, mu_tot, constr);
-
-            const bool do_ene = true;
-            const bool do_log = true;
-            Awh*       awh    = nullptr;
-            const bool do_dr  = ir->nstdisreout != 0;
-            const bool do_or  = ir->nstorireout != 0;
-
-            EnergyOutput::printAnnealingTemperatures(do_log ? fplog : nullptr, groups, &(ir->opts));
-            energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or,
-                                               do_log ? fplog : nullptr, step, t, fr->fcdata.get(), awh);
-
-            if (do_per_step(step, ir->nstlog))
-            {
-                if (fflush(fplog) != 0)
+            if (MASTER(cr))
+            {
+                mimicCommunicator.sendMultipoleOrder(nonConstGlobalTopology);
+            }
+        }
+        else if (command == MCL_RECV_ATOM_COORDS)
+        {
+            wallcycle_start(wcycle, ewcSTEP);
+            if (MASTER(cr))
+            {
+                mimicCommunicator.getCoords(&state_global->x, state_global->natoms);
+            }
+
+            if (ir->efep != efepNO)
+            {
+                state->lambda = currentLambdas(step, *(ir->fepvals), state_global->fep_state);
+            }
+
+            if (MASTER(cr))
+            {
+                const bool constructVsites = ((vsite != nullptr) && mdrunOptions.rerunConstructVsites);
+                if (constructVsites && DOMAINDECOMP(cr))
                 {
-                    gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
+                    gmx_fatal(FARGS,
+                              "Vsite recalculation with -rerun is not implemented with domain "
+                              "decomposition, "
+                              "use a single rank");
+                }
+                if (constructVsites)
+                {
+                    wallcycle_start(wcycle, ewcVSITECONSTR);
+                    vsite->construct(state->x, ir->delta_t, state->v, state->box);
+                    wallcycle_stop(wcycle, ewcVSITECONSTR);
                 }
             }
-        }
 
-        /* Print the remaining wall clock time for the run */
-        if (isMasterSimMasterRank(ms, MASTER(cr)) && (mdrunOptions.verbose || gmx_got_usr_signal()))
+            if (DOMAINDECOMP(cr))
+            {
+                /* Repartition the domain decomposition */
+                const bool bMasterState = true;
+                dd_partition_system(fplog, mdlog, step, cr, bMasterState, nstglobalcomm, state_global,
+                                    *top_global, ir, imdSession, pull_work, state, &f, mdAtoms,
+                                    &top, fr, vsite, constr, nrnb, wcycle, mdrunOptions.verbose);
+                shouldCheckNumberOfBondedInteractions = true;
+            }
+
+            if (MASTER(cr))
+            {
+                EnergyOutput::printHeader(fplog, step, t); /* can we improve the information printed here? */
+            }
+
+            if (ir->efep != efepNO)
+            {
+                update_mdatoms(mdatoms, state->lambda[efptMASS]);
+            }
+        }
+        else if (command == MCL_COMPUTE_FORCE_ENERGY)
         {
+            force_flags = (GMX_FORCE_STATECHANGED | GMX_FORCE_DYNAMICBOX | GMX_FORCE_ALLFORCES
+                           | GMX_FORCE_VIRIAL | // TODO: Get rid of this once #2649 is solved
+                           GMX_FORCE_ENERGY | (doFreeEnergyPerturbation ? GMX_FORCE_DHDL : 0));
+
             if (shellfc)
             {
-                fprintf(stderr, "\n");
+                /* Now is the time to relax the shells */
+                relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, enforcedRotation, step, ir,
+                                    imdSession, pull_work, bNS, force_flags, &top, constr, enerd,
+                                    state->natoms, state->x.arrayRefWithPadding(),
+                                    state->v.arrayRefWithPadding(), state->box, state->lambda,
+                                    &state->hist, &f.view(), force_vir, mdatoms, nrnb, wcycle, shellfc,
+                                    fr, runScheduleWork, t, mu_tot, vsite, ddBalanceRegionHandler);
+            }
+            else
+            {
+                /* The coordinates (x) are shifted (to get whole molecules)
+                 * in do_force.
+                 * This is parallellized as well, and does communication too.
+                 * Check comments in sim_util.c
+                 */
+                Awh*       awh = nullptr;
+                gmx_edsam* ed  = nullptr;
+                do_force(fplog, cr, ms, ir, awh, enforcedRotation, imdSession, pull_work, step, nrnb,
+                         wcycle, &top, state->box, state->x.arrayRefWithPadding(), &state->hist,
+                         &f.view(), force_vir, mdatoms, enerd, state->lambda, fr, runScheduleWork,
+                         vsite, mu_tot, t, ed, GMX_FORCE_NS | force_flags, ddBalanceRegionHandler);
+            }
+
+            /* Now we have the energies and forces corresponding to the
+             * coordinates at time t.
+             */
+            {
+                const bool isCheckpointingStep = false;
+                const bool doRerun             = false;
+                const bool bSumEkinhOld        = false;
+                do_md_trajectory_writing(fplog, cr, nfile, fnm, step, step_rel, t, ir, state,
+                                         state_global, observablesHistory, top_global, fr, outf,
+                                         energyOutput, ekind, f.view().force(), isCheckpointingStep,
+                                         doRerun, isLastStep, mdrunOptions.writeConfout, bSumEkinhOld);
+            }
+
+            stopHandler->setSignal();
+
+            {
+                const bool          doInterSimSignal = false;
+                const bool          doIntraSimSignal = true;
+                bool                bSumEkinhOld     = false;
+                t_vcm*              vcm              = nullptr;
+                SimulationSignaller signaller(&signals, cr, ms, doInterSimSignal, doIntraSimSignal);
+
+                compute_globals(gstat, cr, ir, fr, ekind, makeConstArrayRef(state->x),
+                                makeConstArrayRef(state->v), state->box, mdatoms, nrnb, vcm, wcycle,
+                                enerd, nullptr, nullptr, nullptr, nullptr, constr, &signaller,
+                                state->box, &totalNumberOfBondedInteractions, &bSumEkinhOld,
+                                CGLO_GSTAT | CGLO_ENERGY
+                                        | (shouldCheckNumberOfBondedInteractions ? CGLO_CHECK_NUMBER_OF_BONDED_INTERACTIONS
+                                                                                 : 0));
+                checkNumberOfBondedInteractions(mdlog, cr, totalNumberOfBondedInteractions,
+                                                top_global, &top, makeConstArrayRef(state->x),
+                                                state->box, &shouldCheckNumberOfBondedInteractions);
+            }
+
+            /* Note: this is OK, but there are some numerical precision issues with using the convergence of
+               the virial that should probably be addressed eventually. state->veta has better properies,
+               but what we actually need entering the new cycle is the new shake_vir value. Ideally, we could
+               generate the new shake_vir, but test the veta value for convergence.  This will take some thought. */
+
+            if (ir->efep != efepNO)
+            {
+                /* Sum up the foreign energy and dhdl terms for md and sd.
+                   Currently done every step so that dhdl is correct in the .edr */
+                accumulateKineticLambdaComponents(enerd, state->lambda, *ir->fepvals);
             }
-            print_time(stderr, walltime_accounting, step, ir, cr);
-        }
 
-        cycles = wallcycle_stop(wcycle, ewcSTEP);
-        if (DOMAINDECOMP(cr) && wcycle)
+            /* Output stuff */
+            if (MASTER(cr))
+            {
+                const bool bCalcEnerStep = true;
+                energyOutput.addDataAtEnergyStep(
+                        doFreeEnergyPerturbation, bCalcEnerStep, t, mdatoms->tmass, enerd,
+                        ir->fepvals, ir->expandedvals, state->box,
+                        PTCouplingArrays({ state->boxv, state->nosehoover_xi, state->nosehoover_vxi,
+                                           state->nhpres_xi, state->nhpres_vxi }),
+                        state->fep_state, shake_vir, force_vir, total_vir, pres, ekind, mu_tot, constr);
+
+                const bool do_ene = true;
+                const bool do_log = true;
+                Awh*       awh    = nullptr;
+                const bool do_dr  = ir->nstdisreout != 0;
+                const bool do_or  = ir->nstorireout != 0;
+
+                EnergyOutput::printAnnealingTemperatures(do_log ? fplog : nullptr, groups, &(ir->opts));
+                energyOutput.printStepToEnergyFile(mdoutf_get_fp_ene(outf), do_ene, do_dr, do_or,
+                                                   do_log ? fplog : nullptr, step, t,
+                                                   fr->fcdata.get(), awh);
+
+                if (do_per_step(step, ir->nstlog))
+                {
+                    if (fflush(fplog) != 0)
+                    {
+                        gmx_fatal(FARGS, "Cannot flush logfile - maybe you are out of disk space?");
+                    }
+                }
+            }
+
+            /* Print the remaining wall clock time for the run */
+            if (isMasterSimMasterRank(ms, MASTER(cr)) && (mdrunOptions.verbose || gmx_got_usr_signal()))
+            {
+                if (shellfc)
+                {
+                    fprintf(stderr, "\n");
+                }
+                print_time(stderr, walltime_accounting, step, ir, cr);
+            }
+
+            cycles = wallcycle_stop(wcycle, ewcSTEP);
+            if (DOMAINDECOMP(cr) && wcycle)
+            {
+                dd_cycles_add(cr->dd, cycles, ddCyclStep);
+            }
+
+            /* increase the MD step number */
+            step++;
+            step_rel++;
+        }
+        else
         {
-            dd_cycles_add(cr->dd, cycles, ddCyclStep);
+            gmx_fatal(FARGS, "Unrecognized MiMiC command: %d", command);
         }
-
-        /* increase the MD step number */
-        step++;
-        step_rel++;
     }
     /* End of main MD loop */
 
@@ -573,7 +712,7 @@ void gmx::LegacySimulator::do_mimic()
 
     if (MASTER(cr))
     {
-        MimicCommunicator::finalize();
+        mimicCommunicator.finalize();
     }
 
     if (!thisRankHasDuty(cr, DUTY_PME))
diff --git a/src/gromacs/mdrunutility/threadaffinity.cpp b/src/gromacs/mdrunutility/threadaffinity.cpp
index 8cf377acbf..f299888cd8 100644
--- a/src/gromacs/mdrunutility/threadaffinity.cpp
+++ b/src/gromacs/mdrunutility/threadaffinity.cpp
@@ -542,7 +542,7 @@ static bool detectDefaultAffinityMask(const int nthreads_hw_avail)
     if (mpiIsInitialized)
     {
         bool maskToReduce = detectedDefaultAffinityMask;
-        MPI_Allreduce(&maskToReduce, &detectedDefaultAffinityMask, 1, MPI_C_BOOL, MPI_LAND, MPI_COMM_WORLD);
+        MPI_Allreduce(&maskToReduce, &detectedDefaultAffinityMask, 1, MPI_C_BOOL, MPI_LAND, MPI_WORLD_STUB);
     }
 #endif
 
diff --git a/src/gromacs/mimic/communicator.cpp b/src/gromacs/mimic/communicator.cpp
index 12dcf874d3..28a6087fcf 100644
--- a/src/gromacs/mimic/communicator.cpp
+++ b/src/gromacs/mimic/communicator.cpp
@@ -46,6 +46,7 @@
 #if GMX_MIMIC
 #    include <DataTypes.h>
 #    include <MessageApi.h>
+#    include <Requests.h>
 #endif
 
 // When not built in a configuration with QMMM support, much of this
@@ -102,26 +103,61 @@ void gmx::MimicCommunicator::init()
 {
     char path[GMX_PATH_MAX];
     gmx_getcwd(path, GMX_PATH_MAX);
-    MCL_init_client(path);
+    char del = ';';
+    MCL_handshake(path, &del, 0);
 }
 
-void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVector<gmx::RVec> coords)
+void gmx::MimicCommunicator::sendAtomSpecies(gmx_mtop_t* mtop) {
+    std::vector<int>        atomTypes;
+    atomTypes.reserve(static_cast<size_t>(mtop->natoms));
+
+    int offset = 0;
+    for (const gmx_molblock_t& molblock : mtop->molblock)
+    {
+        gmx_moltype_t* type = &mtop->moltype[molblock.type];
+        for (int mol = 0; mol < molblock.nmol; ++mol)
+        {
+            for (int at = 0; at < type->atoms.nr; ++at)
+            {
+                int  atomtype = type->atoms.atom[at].type;
+                atomTypes.push_back(atomtype + 1);
+            }
+        }
+    }
+    MCL_send(&*atomTypes.begin(), mtop->natoms, TYPE_INT, MCL_DATA, 0);
+}
+void gmx::MimicCommunicator::sendFragmentCount(gmx_mtop_t* mtop)
 {
-    MCL_send(&mtop->natoms, 1, TYPE_INT, 0);
-    MCL_send(&mtop->atomtypes.nr, 1, TYPE_INT, 0);
+    int nMolecules = 0;
+    for (const gmx_molblock_t& molblock : mtop->molblock)
+    {
+        nMolecules += molblock.nmol;
+    }
+    MCL_send(&nMolecules, 1, TYPE_INT, MCL_DATA, 0);
+}
 
-    std::vector<int>        atomTypes;
-    std::vector<int>        nAtomsMol;
-    std::vector<int>        idOrder;
-    std::vector<double>     charges;
-    std::vector<double>     masses(mtop->atomtypes.nr, -1);
-    std::vector<int>        elements(mtop->atomtypes.nr, -1);
-    std::vector<int>        bonds;
-    std::vector<double>     bondLengths;
-    std::unordered_set<int> existingTypes;
+void gmx::MimicCommunicator::sendBondsNumber(gmx_mtop_t *mtop)
+{
+    int bondsNumber = 0;
+    for (const gmx_molblock_t& molblock : mtop->molblock)
+    {
+        gmx_moltype_t* type = &mtop->moltype[molblock.type];
+        for (int mol = 0; mol < molblock.nmol; ++mol)
+        {
+            int nconstr  = type->ilist[F_CONSTR].size() / 3;
+            int nconstrc = type->ilist[F_CONSTRNC].size() / 3;
+            int nsettle  = type->ilist[F_SETTLE].size() / 4;
 
-    atomTypes.reserve(static_cast<size_t>(mtop->natoms));
-    charges.reserve(static_cast<size_t>(mtop->natoms));
+            bondsNumber += nconstr + nconstrc + nsettle * 3;
+        }
+    }
+
+    MCL_send(&bondsNumber, 1, TYPE_INT, MCL_DATA, 0);
+}
+
+void gmx::MimicCommunicator::sendBondAtoms(gmx_mtop_t *mtop)
+{
+    std::vector<int>        bonds;
 
     int offset = 0;
     for (const gmx_molblock_t& molblock : mtop->molblock)
@@ -135,13 +171,10 @@ void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVector<gmx
 
             for (int ncon = 0; ncon < nconstr + nconstrc; ++ncon)
             {
-                int contype = type->ilist[F_CONSTR].iatoms[0];
                 int at1     = type->ilist[F_CONSTR].iatoms[1];
                 int at2     = type->ilist[F_CONSTR].iatoms[2];
                 bonds.push_back(offset + at1 + 1);
                 bonds.push_back(offset + at2 + 1);
-                bondLengths.push_back(static_cast<double>(mtop->ffparams.iparams[contype].constr.dA)
-                                      / BOHR2NM);
             }
 
             for (int ncon = 0; ncon < nsettle; ++ncon)
@@ -150,8 +183,6 @@ void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVector<gmx
                 t_iatom h1;
                 t_iatom h2;
 
-                int contype = type->ilist[F_SETTLE].iatoms[0];
-
                 ox = type->ilist[F_SETTLE].iatoms[1];
                 h1 = type->ilist[F_SETTLE].iatoms[2];
                 h2 = type->ilist[F_SETTLE].iatoms[3];
@@ -164,6 +195,38 @@ void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVector<gmx
 
                 bonds.push_back(offset + h1 + 1);
                 bonds.push_back(offset + h2 + 1);
+            }
+
+            offset += type->atoms.nr;
+        }
+    }
+
+    MCL_send(&*bonds.begin(), bonds.size(), TYPE_INT, MCL_DATA, 0);
+}
+
+void gmx::MimicCommunicator::sendBondLengths(gmx_mtop_t* mtop)
+{
+    std::vector<double>     bondLengths;
+
+    for (const gmx_molblock_t& molblock : mtop->molblock)
+    {
+        gmx_moltype_t* type = &mtop->moltype[molblock.type];
+        for (int mol = 0; mol < molblock.nmol; ++mol)
+        {
+            int nconstr  = type->ilist[F_CONSTR].size() / 3;
+            int nconstrc = type->ilist[F_CONSTRNC].size() / 3;
+            int nsettle  = type->ilist[F_SETTLE].size() / 4;
+
+            for (int ncon = 0; ncon < nconstr + nconstrc; ++ncon)
+            {
+                int contype = type->ilist[F_CONSTR].iatoms[0];
+                bondLengths.push_back(static_cast<double>(mtop->ffparams.iparams[contype].constr.dA)
+                                      / BOHR2NM);
+            }
+
+            for (int ncon = 0; ncon < nsettle; ++ncon)
+            {
+                int contype = type->ilist[F_SETTLE].iatoms[0];
                 bondLengths.push_back(static_cast<double>(mtop->ffparams.iparams[contype].constr.dA)
                                       / BOHR2NM);
                 bondLengths.push_back(static_cast<double>(mtop->ffparams.iparams[contype].constr.dA)
@@ -172,65 +235,128 @@ void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVector<gmx
                                       / BOHR2NM);
             }
 
-            nAtomsMol.push_back(type->atoms.nr);
+        }
+    }
+    MCL_send(&*bondLengths.begin(), bondLengths.size(), TYPE_DOUBLE, MCL_DATA, 0);
+}
+
+void gmx::MimicCommunicator::sendAnglesNumber(gmx_mtop_t* mtop)
+{
+    int dummyAngles = 0;
+    MCL_send(&dummyAngles, 1, TYPE_INT, MCL_DATA, 0);
+}
+
+void gmx::MimicCommunicator::sendAtomElements(gmx_mtop_t* mtop)
+{
+    std::vector<int>        elements(mtop->atomtypes.nr, -1);
+    std::unordered_set<int> existingTypes;
+
+    for (const gmx_molblock_t& molblock : mtop->molblock)
+    {
+        gmx_moltype_t* type = &mtop->moltype[molblock.type];
+        for (int mol = 0; mol < molblock.nmol; ++mol)
+        {
             for (int at = 0; at < type->atoms.nr; ++at)
             {
                 int  atomtype = type->atoms.atom[at].type;
-                auto charge   = static_cast<double>(type->atoms.atom[at].q);
-                idOrder.push_back(offset + 1);
-                offset++;
-                atomTypes.push_back(atomtype + 1);
-                charges.push_back(charge);
                 if (existingTypes.insert(atomtype).second)
                 {
-                    masses[atomtype]   = type->atoms.atom[at].m;
                     elements[atomtype] = type->atoms.atom[at].atomnumber;
                 }
             }
         }
     }
-    // sending atom types
-    MCL_send(&*atomTypes.begin(), mtop->natoms, TYPE_INT, 0);
 
-    int max_multipole_order = 0;
-    // sending multipole orders
-    MCL_send(&max_multipole_order, 1, TYPE_INT, 0);
-
-    int nMolecules = nAtomsMol.size();
-    // sending molecule number
-    MCL_send(&nMolecules, 1, TYPE_INT, 0);
+    MCL_send(&*elements.begin(), mtop->atomtypes.nr, TYPE_INT, MCL_DATA, 0);
+}
 
-    // sending number of atoms per molecules
-    MCL_send(&*nAtomsMol.begin(), nAtomsMol.size(), TYPE_INT, 0);
+void gmx::MimicCommunicator::sendAtomMasses(gmx_mtop_t* mtop)
+{
+    std::vector<double>     masses(mtop->atomtypes.nr, -1);
+    std::unordered_set<int> existingTypes;
 
-    int nBonds = bonds.size() / 2;
-    // sending number of bond constraints
-    MCL_send(&nBonds, 1, TYPE_INT, 0);
+    for (const gmx_molblock_t& molblock : mtop->molblock)
+    {
+        gmx_moltype_t* type = &mtop->moltype[molblock.type];
+        for (int mol = 0; mol < molblock.nmol; ++mol)
+        {
+            for (int at = 0; at < type->atoms.nr; ++at)
+            {
+                int  atomtype = type->atoms.atom[at].type;
+                if (existingTypes.insert(atomtype).second)
+                {
+                    masses[atomtype]   = type->atoms.atom[at].m;
+                }
+            }
+        }
+    }
 
-    // sending number of angle constraints
-    MCL_send(&max_multipole_order, 1, TYPE_INT, 0);
+    MCL_send(&*masses.begin(), masses.size(), TYPE_DOUBLE, MCL_DATA, 0);
+}
 
-    if (nBonds > 0)
+void gmx::MimicCommunicator::sendFragmentAtomNum(gmx_mtop_t* mtop)
+{
+    std::vector<int>        nAtomsMol;
+    for (const gmx_molblock_t& molblock : mtop->molblock)
     {
-        // sending bonded atoms indices
-        MCL_send(&*bonds.begin(), bonds.size(), TYPE_INT, 0);
-
-        // sending bond lengths
-        MCL_send(&*bondLengths.begin(), bondLengths.size(), TYPE_DOUBLE, 0);
+        gmx_moltype_t* type = &mtop->moltype[molblock.type];
+        for (int mol = 0; mol < molblock.nmol; ++mol)
+        {
+            nAtomsMol.push_back(type->atoms.nr);
+        }
     }
+    // sending number of atoms per molecules
+    MCL_send(&*nAtomsMol.begin(), nAtomsMol.size(), TYPE_INT, MCL_DATA, 0);
+}
+
+void gmx::MimicCommunicator::sendFragmentAtomIds(gmx_mtop_t* mtop)
+{
+    std::vector<int>        idOrder;
 
-    // sending array of atomic charges
-    MCL_send(&*charges.begin(), mtop->natoms, TYPE_DOUBLE, 0);
+    int offset = 0;
+    for (const gmx_molblock_t& molblock : mtop->molblock)
+    {
+        gmx_moltype_t* type = &mtop->moltype[molblock.type];
+        for (int mol = 0; mol < molblock.nmol; ++mol)
+        {
+            for (int at = 0; at < type->atoms.nr; ++at)
+            {
+                idOrder.push_back(offset + 1);
+                offset++;
+            }
+        }
+    }
+    MCL_send(&*idOrder.begin(), idOrder.size(), TYPE_INT, MCL_DATA, 0);
+}
 
-    // sending array of atomic masses
-    MCL_send(&*masses.begin(), mtop->atomtypes.nr, TYPE_DOUBLE, 0);
+void gmx::MimicCommunicator::sendMultipoleOrder(gmx_mtop_t* mtop)
+{
+    int multipleOrder = 0;
+    MCL_send(&multipleOrder, 1, TYPE_INT, MCL_DATA, 0);
+}
 
-    // sending ids of atoms per molecule
-    MCL_send(&*idOrder.begin(), idOrder.size(), TYPE_INT, 0);
+void gmx::MimicCommunicator::sendCharges(gmx_mtop_t* mtop)
+{
+    std::vector<double>     charges;
+    charges.reserve(static_cast<size_t>(mtop->natoms));
 
-    // sending list of elements
-    MCL_send(&*elements.begin(), mtop->atomtypes.nr, TYPE_INT, 0);
+    for (const gmx_molblock_t& molblock : mtop->molblock)
+    {
+        gmx_moltype_t* type = &mtop->moltype[molblock.type];
+        for (int mol = 0; mol < molblock.nmol; ++mol)
+        {
+            for (int at = 0; at < type->atoms.nr; ++at)
+            {
+                auto charge   = static_cast<double>(type->atoms.atom[at].q);
+                charges.push_back(charge);
+            }
+        }
+    }
+    MCL_send(&*charges.begin(), mtop->natoms, TYPE_DOUBLE, MCL_DATA, 0);
+}
 
+void gmx::MimicCommunicator::sendCoords(PaddedHostVector<gmx::RVec> coords)
+{
     std::vector<double> convertedCoords;
     for (auto& coord : coords)
     {
@@ -240,20 +366,27 @@ void gmx::MimicCommunicator::sendInitData(gmx_mtop_t* mtop, PaddedHostVector<gmx
     }
 
     // sending array of coordinates
-    MCL_send(&*convertedCoords.begin(), 3 * mtop->natoms, TYPE_DOUBLE, 0);
+    MCL_send(&*convertedCoords.begin(), convertedCoords.size(), TYPE_DOUBLE, MCL_DATA, 0);
 }
 
-int64_t gmx::MimicCommunicator::getStepNumber()
-{
-    int steps;
-    MCL_receive(&steps, 1, TYPE_INT, 0);
-    return steps;
+void gmx::MimicCommunicator::sendAtomNum(gmx_mtop_t* mtop) {
+    MCL_send(&mtop->natoms, 1, TYPE_INT, MCL_DATA, 0);
+}
+
+void gmx::MimicCommunicator::sendSpeciesNum(gmx_mtop_t* mtop) {
+    MCL_send(&mtop->atomtypes.nr, 1, TYPE_INT, MCL_DATA, 0);
+}
+
+int gmx::MimicCommunicator::getCommand() {
+    int command{};
+    MCL_receive(&command, 1, TYPE_INT, MCL_COMMAND, 0);
+    return command;
 }
 
 void gmx::MimicCommunicator::getCoords(PaddedHostVector<RVec>* x, const int natoms)
 {
     std::vector<double> coords(natoms * 3);
-    MCL_receive(&*coords.begin(), 3 * natoms, TYPE_DOUBLE, 0);
+    MCL_receive(&*coords.begin(), 3 * natoms, TYPE_DOUBLE, MCL_DATA, 0);
     for (int j = 0; j < natoms; ++j)
     {
         (*x)[j][0] = static_cast<real>(coords[j * 3] * BOHR2NM);
@@ -265,7 +398,7 @@ void gmx::MimicCommunicator::getCoords(PaddedHostVector<RVec>* x, const int nato
 void gmx::MimicCommunicator::sendEnergies(real energy)
 {
     double convertedEnergy = energy / (HARTREE2KJ * AVOGADRO);
-    MCL_send(&convertedEnergy, 1, TYPE_DOUBLE, 0);
+    MCL_send(&convertedEnergy, 1, TYPE_DOUBLE, MCL_DATA, 0);
 }
 
 void gmx::MimicCommunicator::sendForces(gmx::ArrayRef<gmx::RVec> forces, int natoms)
@@ -277,7 +410,7 @@ void gmx::MimicCommunicator::sendForces(gmx::ArrayRef<gmx::RVec> forces, int nat
         convertedForce.push_back(static_cast<real>(forces[j][1]) / HARTREE_BOHR2MD);
         convertedForce.push_back(static_cast<real>(forces[j][2]) / HARTREE_BOHR2MD);
     }
-    MCL_send(&*convertedForce.begin(), convertedForce.size(), TYPE_DOUBLE, 0);
+    MCL_send(&*convertedForce.begin(), convertedForce.size(), TYPE_DOUBLE, MCL_DATA, 0);
 }
 
 void gmx::MimicCommunicator::finalize()
diff --git a/src/gromacs/mimic/communicator.h b/src/gromacs/mimic/communicator.h
index 6d2b9c2761..436a13c784 100644
--- a/src/gromacs/mimic/communicator.h
+++ b/src/gromacs/mimic/communicator.h
@@ -60,27 +60,7 @@ public:
     /*! \brief
      * Initializes the communicator
      */
-    static void init();
-
-    /*! \brief
-     * Sends the data needed for MiMiC initialization
-     *
-     * That includes number of atoms, element numbers, charges, masses,
-     * maximal order of multipoles (0 for point-charge forcefields),
-     * number of molecules, number of atoms per each molecule,
-     * bond constraints data
-     *
-     * @param mtop global topology data
-     * @param coords coordinates of all atoms
-     */
-    static void sendInitData(gmx_mtop_t* mtop, PaddedHostVector<gmx::RVec> coords);
-
-    /*! \brief
-     * Gets the number of MD steps to perform from MiMiC
-     *
-     * @return nsteps the number of MD steps to perform
-     */
-    static int64_t getStepNumber();
+    void init();
 
     /*! \brief
      * Receive and array of updated atomic coordinates from MiMiC
@@ -88,14 +68,14 @@ public:
      * @param x array of coordinates to fill
      * @param natoms number of atoms in the system
      */
-    static void getCoords(PaddedHostVector<RVec>* x, int natoms);
+    void getCoords(PaddedHostVector<RVec>* x, int natoms);
 
     /*! \brief
      * Send the potential energy value to MiMiC
      *
      * @param energy energy value to send
      */
-    static void sendEnergies(real energy);
+    void sendEnergies(real energy);
 
     /*! \brief
      * Send classical forces acting on all atoms in the system
@@ -104,12 +84,44 @@ public:
      * @param forces array of forces to send
      * @param natoms number of atoms in the system
      */
-    static void sendForces(ArrayRef<gmx::RVec> forces, int natoms);
+    void sendForces(ArrayRef<gmx::RVec> forces, int natoms);
 
     /*! \brief
      * Finish communications and disconnect from the server
      */
-    static void finalize();
+    void finalize();
+
+    void sendAtomSpecies(gmx_mtop_t *mtop);
+
+    void sendFragmentCount(gmx_mtop_t *mtop);
+
+    void sendBondsNumber(gmx_mtop_t *mtop);
+
+    void sendBondAtoms(gmx_mtop_t *mtop);
+
+    void sendBondLengths(gmx_mtop_t *mtop);
+
+    void sendAnglesNumber(gmx_mtop_t *mtop);
+
+    void sendAtomElements(gmx_mtop_t *mtop);
+
+    void sendAtomMasses(gmx_mtop_t *mtop);
+
+    void sendFragmentAtomNum(gmx_mtop_t *mtop);
+
+    void sendFragmentAtomIds(gmx_mtop_t *mtop);
+
+    void sendMultipoleOrder(gmx_mtop_t *mtop);
+
+    void sendCharges(gmx_mtop_t *mtop);
+
+    void sendCoords(PaddedHostVector <gmx::RVec> coords);
+
+    int getCommand();
+
+    void sendAtomNum(gmx_mtop_t *mtop);
+
+    void sendSpeciesNum(gmx_mtop_t *mtop);
 };
 
 } // namespace gmx
diff --git a/src/gromacs/tools/pme_error.cpp b/src/gromacs/tools/pme_error.cpp
index ba1fb3b142..7f2c31fecc 100644
--- a/src/gromacs/tools/pme_error.cpp
+++ b/src/gromacs/tools/pme_error.cpp
@@ -1121,7 +1121,7 @@ int gmx_pme_error(int argc, char* argv[])
 
 #define NFILE asize(fnm)
 
-    CommrecHandle commrecHandle = init_commrec(MPI_COMM_WORLD);
+    CommrecHandle commrecHandle = init_commrec(MPI_WORLD_STUB);
     t_commrec*    cr            = commrecHandle.get();
     PCA_Flags                   = PCA_NOEXIT_ON_ARGS;
 
diff --git a/src/gromacs/utility/basenetwork.cpp b/src/gromacs/utility/basenetwork.cpp
index 6fade03e57..bae56d5259 100644
--- a/src/gromacs/utility/basenetwork.cpp
+++ b/src/gromacs/utility/basenetwork.cpp
@@ -74,7 +74,7 @@ int gmx_node_num()
     }
 #    endif
     int i;
-    (void)MPI_Comm_size(MPI_COMM_WORLD, &i);
+    (void)MPI_Comm_size(MPI_WORLD_STUB, &i);
     return i;
 #endif
 }
@@ -91,7 +91,7 @@ int gmx_node_rank()
     }
 #    endif
     int i;
-    (void)MPI_Comm_rank(MPI_COMM_WORLD, &i);
+    (void)MPI_Comm_rank(MPI_WORLD_STUB, &i);
     return i;
 #endif
 }
@@ -153,7 +153,7 @@ int gmx_physicalnode_id_hash()
 void gmx_broadcast_world(int size, void* buffer)
 {
 #if GMX_MPI
-    MPI_Bcast(buffer, size, MPI_BYTE, 0, MPI_COMM_WORLD);
+    MPI_Bcast(buffer, size, MPI_BYTE, 0, MPI_WORLD_STUB);
 #else
     GMX_UNUSED_VALUE(size);
     GMX_UNUSED_VALUE(buffer);
@@ -163,7 +163,7 @@ void gmx_broadcast_world(int size, void* buffer)
 #if GMX_LIB_MPI
 void gmx_abort(int errorno)
 {
-    MPI_Abort(MPI_COMM_WORLD, errorno);
+    MPI_Abort(MPI_WORLD_STUB, errorno);
     std::abort();
 }
 #endif
diff --git a/src/gromacs/utility/fatalerror.cpp b/src/gromacs/utility/fatalerror.cpp
index f29b400a80..d09225c000 100644
--- a/src/gromacs/utility/fatalerror.cpp
+++ b/src/gromacs/utility/fatalerror.cpp
@@ -186,7 +186,7 @@ void gmx_exit_on_fatal_error(ExitType exitType, int returnValue)
             case ExitType_NonMasterAbort:
                 // Let all other processes wait till the master has printed
                 // the error message and issued MPI_Abort.
-                MPI_Barrier(MPI_COMM_WORLD);
+                MPI_Barrier(MPI_WORLD_STUB);
                 break;
         }
     }
diff --git a/src/gromacs/utility/gmxmpi.h b/src/gromacs/utility/gmxmpi.h
index 5c2a247526..f226e3644c 100644
--- a/src/gromacs/utility/gmxmpi.h
+++ b/src/gromacs/utility/gmxmpi.h
@@ -60,6 +60,7 @@
 /* disable bindings for SGI MPT also */
 #    define MPI_NO_CPPBIND 1
 #    include <mpi.h>
+extern MPI_Comm MPI_WORLD_STUB;
 /* Starting with 2.2 MPI_INT64_T is required. Earlier version still might have it.
    In theory MPI_Datatype doesn't have to be a #define, but current available MPI
    implementations (OpenMPI + MPICH (+derivates)) use #define and future versions
diff --git a/src/gromacs/utility/init.cpp b/src/gromacs/utility/init.cpp
index 1f6e07414b..f2ca2a3986 100644
--- a/src/gromacs/utility/init.cpp
+++ b/src/gromacs/utility/init.cpp
@@ -53,6 +53,10 @@
 #    include "gromacs/utility/gmxmpi.h"
 #endif
 
+#include "MessageApi.h"
+
+MPI_Comm MPI_WORLD_STUB = MPI_COMM_WORLD;
+
 namespace gmx
 {
 
@@ -114,6 +118,7 @@ void init(int* argc, char*** argv) // NOLINT(readability-non-const-parameter)
 #    else
         MPI_Init(argc, argv);
 #    endif
+        MCL_init(&MPI_WORLD_STUB);
     }
     // Bump the counter to record this initialization event
     g_initializationCounter++;
@@ -137,7 +142,7 @@ void finalize()
          * with buggy MPI implementations that could cause
          * unfinished processes to terminate.
          */
-        MPI_Barrier(MPI_COMM_WORLD);
+        MPI_Barrier(MPI_WORLD_STUB);
 
         /* Apparently certain mpich implementations cause problems
          * with MPI_Finalize. In that case comment out MPI_Finalize.
diff --git a/src/programs/mdrun/mdrun.cpp b/src/programs/mdrun/mdrun.cpp
index 77a25f1fdf..73ded62fd4 100644
--- a/src/programs/mdrun/mdrun.cpp
+++ b/src/programs/mdrun/mdrun.cpp
@@ -70,6 +70,10 @@
 #include "gromacs/utility/basenetwork.h"
 #include "gromacs/utility/physicalnodecommunicator.h"
 
+#if GMX_MIMIC
+#    include <MessageApi.h>
+#endif
+
 #include "mdrun_main.h"
 
 namespace gmx
@@ -79,7 +83,7 @@ int gmx_mdrun(int argc, char* argv[])
 {
     // Set up the communicator, where possible (see docs for
     // SimulationContext).
-    MPI_Comm                 communicator = GMX_LIB_MPI ? MPI_COMM_WORLD : MPI_COMM_NULL;
+    MPI_Comm                 communicator = GMX_LIB_MPI ? MPI_WORLD_STUB : MPI_COMM_NULL;
     PhysicalNodeCommunicator physicalNodeCommunicator(communicator, gmx_physicalnode_id_hash());
     std::unique_ptr<gmx_hw_info_t> hwinfo = gmx_detect_hardware(physicalNodeCommunicator);
     return gmx_mdrun(communicator, *hwinfo, argc, argv);
