Merge branch 'release-2018'
authorAleksei Iupinov <a.yupinov@gmail.com>
Tue, 27 Feb 2018 13:08:17 +0000 (14:08 +0100)
committerAleksei Iupinov <a.yupinov@gmail.com>
Tue, 27 Feb 2018 14:05:53 +0000 (15:05 +0100)
Conflict: multisim barrier handling in src/programs/mdrun/runner.cpp

Amended the test introduced in I16a3705a0d59a to exclude the implicit
solvent output which was disabled in Ib241555ff3d8e60.

Change-Id: Idf382c9c7ba92439097fcc9008cb9cc813daaab7

1  2 
docs/CMakeLists.txt
src/gromacs/gmxpreprocess/tests/readir.cpp
src/gromacs/gmxpreprocess/tests/refdata/GetIrTest_DefineHandlesAssignmentOnRhs.xml
src/gromacs/mdlib/update.cpp
src/programs/mdrun/runner.cpp

diff --combined docs/CMakeLists.txt
index e8114880a0fc60081dc7c8b081912d59619ed007,67193bf553a0586746fba1c9df9f4514719a5b9b..f18d56fd4a1cab6158aaef5a8ccb555d9b7e7135
@@@ -106,7 -106,6 +106,7 @@@ if (SPHINX_FOUND
          dev-manual/build-system.rst
          dev-manual/commitstyle.rst
          dev-manual/documentation-generation.rst
 +        dev-manual/contribute.rst
          dev-manual/doxygen.rst
          dev-manual/error-handling.rst
          dev-manual/formatting.rst
          release-notes/2018/major/removed-features.rst
          release-notes/2018/major/portability.rst
          release-notes/2018/major/miscellaneous.rst
+         release-notes/2016/2016.5.rst
+         release-notes/2016/2016.4.rst
+         release-notes/2016/2016.3.rst
+         release-notes/2016/2016.2.rst
+         release-notes/2016/2016.1.rst
+         release-notes/2016/major/highlights.rst
+         release-notes/2016/major/new-features.rst
+         release-notes/2016/major/performance.rst
+         release-notes/2016/major/tools.rst
+         release-notes/2016/major/bugs-fixed.rst
+         release-notes/2016/major/removed-features.rst
+         release-notes/2016/major/miscellaneous.rst
+         release-notes/older/index.rst
          user-guide/index.rst
          user-guide/cutoff-schemes.rst
          user-guide/getting-started.rst
index 4b72470dec36f67b769cde621b5124d91d018a37,73768107ee8a46ad7e9757ae4eedb55fb58d53f3..6c07d2eab4ea9ae010d55de1195806bb1d15a015
@@@ -71,7 -71,8 +71,7 @@@ namespace tes
  class GetIrTest : public ::testing::Test
  {
      public:
 -        GetIrTest() : fileManager_(), data_(), checker_(data_.rootChecker()),
 -                      ir_(), mdModules_(), opts_(),
 +        GetIrTest() : fileManager_(), ir_(), mdModules_(), opts_(),
                        wi_(init_warning(FALSE, 0)), wiGuard_(wi_)
  
          {
  
              get_ir(inputMdpFilename.c_str(), outputMdpFilename.c_str(),
                     &mdModules_, &ir_, &opts_, WriteMdpHeader::no, wi_);
 -            bool failure = warning_errors_exist(wi_);
 -            checker_.checkBoolean(failure, "Error parsing mdp file");
 +
 +            // Now check
 +            bool                 failure = warning_errors_exist(wi_);
 +            TestReferenceData    data;
 +            TestReferenceChecker checker(data.rootChecker());
 +            checker.checkBoolean(failure, "Error parsing mdp file");
              warning_reset(wi_);
  
              auto outputMdpContents = TextReader::readFileToString(outputMdpFilename);
 -            checker_.checkString(outputMdpContents, "OutputMdpFile");
 +            checker.checkString(outputMdpContents, "OutputMdpFile");
          }
  
          TestFileManager                    fileManager_;
 -        TestReferenceData                  data_;
 -        TestReferenceChecker               checker_;
          t_inputrec                         ir_;
          MDModules                          mdModules_;
          t_gromppopts                       opts_;
@@@ -157,6 -156,14 +157,14 @@@ TEST_F(GetIrTest, UserErrorsSilentlyTol
      runTest(joinStrings(inputMdpFile, "\n"));
  }
  
+ TEST_F(GetIrTest, DefineHandlesAssignmentOnRhs)
+ {
+     const char *inputMdpFile[] = {
+         "define = -DBOOL -DVAR=VALUE",
+     };
+     runTest(joinStrings(inputMdpFile, "\n"));
+ }
  TEST_F(GetIrTest, EmptyInputWorks)
  {
      const char *inputMdpFile = "";
@@@ -183,23 -190,5 +191,23 @@@ TEST_F(GetIrTest, ProducesOutputFromEle
      runTest(inputMdpFile);
  }
  
 +TEST_F(GetIrTest, TerminatesOnDuplicateOldAndNewKeys)
 +{
 +    const char *inputMdpFile[] = {"verlet-buffer-drift = 1.3", "verlet-buffer-tolerance = 2.7"};
 +    EXPECT_DEATH_IF_SUPPORTED(runTest(joinStrings(inputMdpFile, "\n")), "A parameter is present with both");
 +}
 +
 +TEST_F(GetIrTest, ImplicitSolventNoWorks)
 +{
 +    const char *inputMdpFile = "implicit-solvent = no";
 +    runTest(inputMdpFile);
 +}
 +
 +TEST_F(GetIrTest, ImplicitSolventYesWorks)
 +{
 +    const char *inputMdpFile = "implicit-solvent = yes";
 +    EXPECT_DEATH_IF_SUPPORTED(runTest(inputMdpFile), "Invalid enum");
 +}
 +
  } // namespace
  } // namespace
index 0000000000000000000000000000000000000000,7e032a588acbab1f2d2f4354e84835e654ba7777..95038d7d536f3e7ec5f13b1838b7add1f9dc0ebb
mode 000000,100644..100644
--- /dev/null
@@@ -1,0 -1,345 +1,321 @@@
 -
 -; IMPLICIT SOLVENT ALGORITHM
 -implicit-solvent         = No
 -
 -; GENERALIZED BORN ELECTROSTATICS
 -; Algorithm for calculating Born radii
 -gb-algorithm             = Still
 -; Frequency of calculating the Born radii inside rlist
 -nstgbradii               = 1
 -; Cutoff for Born radii calculation; the contribution from atoms
 -; between rlist and rgbradii is updated every nstlist steps
 -rgbradii                 = 1
 -; Dielectric coefficient of the implicit solvent
 -gb-epsilon-solvent       = 80
 -; Salt concentration in M for Generalized Born models
 -gb-saltconc              = 0
 -; Scaling factors used in the OBC GB model. Default values are OBC(II)
 -gb-obc-alpha             = 1
 -gb-obc-beta              = 0.8
 -gb-obc-gamma             = 4.85
 -gb-dielectric-offset     = 0.009
 -sa-algorithm             = Ace-approximation
 -; Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA
 -; The value -1 will set default value for Still/HCT/OBC GB-models.
 -sa-surface-tension       = -1
+ <?xml version="1.0"?>
+ <?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+ <ReferenceData>
+   <Bool Name="Error parsing mdp file">false</Bool>
+   <String Name="OutputMdpFile">
+ ; VARIOUS PREPROCESSING OPTIONS
+ ; Preprocessor information: use cpp syntax.
+ ; e.g.: -I/home/joe/doe -I/home/mary/roe
+ include                  = 
+ ; e.g.: -DPOSRES -DFLEXIBLE (note these variable names are case sensitive)
+ define                   = -DBOOL -DVAR=VALUE
+ ; RUN CONTROL PARAMETERS
+ integrator               = md
+ ; Start time and timestep in ps
+ tinit                    = 0
+ dt                       = 0.001
+ nsteps                   = 0
+ ; For exact run continuation or redoing part of a run
+ init-step                = 0
+ ; Part index is updated automatically on checkpointing (keeps files separate)
+ simulation-part          = 1
+ ; mode for center of mass motion removal
+ comm-mode                = Linear
+ ; number of steps for center of mass motion removal
+ nstcomm                  = 100
+ ; group(s) for center of mass motion removal
+ comm-grps                = 
+ ; LANGEVIN DYNAMICS OPTIONS
+ ; Friction coefficient (amu/ps) and random seed
+ bd-fric                  = 0
+ ld-seed                  = -1
+ ; ENERGY MINIMIZATION OPTIONS
+ ; Force tolerance and initial step-size
+ emtol                    = 10
+ emstep                   = 0.01
+ ; Max number of iterations in relax-shells
+ niter                    = 20
+ ; Step size (ps^2) for minimization of flexible constraints
+ fcstep                   = 0
+ ; Frequency of steepest descents steps when doing CG
+ nstcgsteep               = 1000
+ nbfgscorr                = 10
+ ; TEST PARTICLE INSERTION OPTIONS
+ rtpi                     = 0.05
+ ; OUTPUT CONTROL OPTIONS
+ ; Output frequency for coords (x), velocities (v) and forces (f)
+ nstxout                  = 0
+ nstvout                  = 0
+ nstfout                  = 0
+ ; Output frequency for energies to log file and energy file
+ nstlog                   = 1000
+ nstcalcenergy            = 100
+ nstenergy                = 1000
+ ; Output frequency and precision for .xtc file
+ nstxout-compressed       = 0
+ compressed-x-precision   = 1000
+ ; This selects the subset of atoms for the compressed
+ ; trajectory file. You can select multiple groups. By
+ ; default, all atoms will be written.
+ compressed-x-grps        = 
+ ; Selection of energy groups
+ energygrps               = 
+ ; NEIGHBORSEARCHING PARAMETERS
+ ; cut-off scheme (Verlet: particle based cut-offs, group: using charge groups)
+ cutoff-scheme            = Verlet
+ ; nblist update frequency
+ nstlist                  = 10
+ ; ns algorithm (simple or grid)
+ ns-type                  = Grid
+ ; Periodic boundary conditions: xyz, no, xy
+ pbc                      = xyz
+ periodic-molecules       = no
+ ; Allowed energy error due to the Verlet buffer in kJ/mol/ps per atom,
+ ; a value of -1 means: use rlist
+ verlet-buffer-tolerance  = 0.005
+ ; nblist cut-off        
+ rlist                    = 1
+ ; long-range cut-off for switched potentials
+ ; OPTIONS FOR ELECTROSTATICS AND VDW
+ ; Method for doing electrostatics
+ coulombtype              = Cut-off
+ coulomb-modifier         = Potential-shift-Verlet
+ rcoulomb-switch          = 0
+ rcoulomb                 = 1
+ ; Relative dielectric constant for the medium and the reaction field
+ epsilon-r                = 1
+ epsilon-rf               = 0
+ ; Method for doing Van der Waals
+ vdw-type                 = Cut-off
+ vdw-modifier             = Potential-shift-Verlet
+ ; cut-off lengths       
+ rvdw-switch              = 0
+ rvdw                     = 1
+ ; Apply long range dispersion corrections for Energy and Pressure
+ DispCorr                 = No
+ ; Extension of the potential lookup tables beyond the cut-off
+ table-extension          = 1
+ ; Separate tables between energy group pairs
+ energygrp-table          = 
+ ; Spacing for the PME/PPPM FFT grid
+ fourierspacing           = 0.12
+ ; FFT grid size, when a value is 0 fourierspacing will be used
+ fourier-nx               = 0
+ fourier-ny               = 0
+ fourier-nz               = 0
+ ; EWALD/PME/PPPM parameters
+ pme-order                = 4
+ ewald-rtol               = 1e-05
+ ewald-rtol-lj            = 0.001
+ lj-pme-comb-rule         = Geometric
+ ewald-geometry           = 3d
+ epsilon-surface          = 0
++implicit-solvent         = no
+ ; OPTIONS FOR WEAK COUPLING ALGORITHMS
+ ; Temperature coupling  
+ tcoupl                   = No
+ nsttcouple               = -1
+ nh-chain-length          = 10
+ print-nose-hoover-chain-variables = no
+ ; Groups to couple separately
+ tc-grps                  = 
+ ; Time constant (ps) and reference temperature (K)
+ tau-t                    = 
+ ref-t                    = 
+ ; pressure coupling     
+ pcoupl                   = No
+ pcoupltype               = Isotropic
+ nstpcouple               = -1
+ ; Time constant (ps), compressibility (1/bar) and reference P (bar)
+ tau-p                    = 1
+ compressibility          = 
+ ref-p                    = 
+ ; Scaling of reference coordinates, No, All or COM
+ refcoord-scaling         = No
+ ; OPTIONS FOR QMMM calculations
+ QMMM                     = no
+ ; Groups treated Quantum Mechanically
+ QMMM-grps                = 
+ ; QM method             
+ QMmethod                 = 
+ ; QMMM scheme           
+ QMMMscheme               = normal
+ ; QM basisset           
+ QMbasis                  = 
+ ; QM charge             
+ QMcharge                 = 
+ ; QM multiplicity       
+ QMmult                   = 
+ ; Surface Hopping       
+ SH                       = 
+ ; CAS space options     
+ CASorbitals              = 
+ CASelectrons             = 
+ SAon                     = 
+ SAoff                    = 
+ SAsteps                  = 
+ ; Scale factor for MM charges
+ MMChargeScaleFactor      = 1
+ ; SIMULATED ANNEALING  
+ ; Type of annealing for each temperature group (no/single/periodic)
+ annealing                = 
+ ; Number of time points to use for specifying annealing in each group
+ annealing-npoints        = 
+ ; List of times at the annealing points for each group
+ annealing-time           = 
+ ; Temp. at each annealing point, for each group.
+ annealing-temp           = 
+ ; GENERATE VELOCITIES FOR STARTUP RUN
+ gen-vel                  = no
+ gen-temp                 = 300
+ gen-seed                 = -1
+ ; OPTIONS FOR BONDS    
+ constraints              = none
+ ; Type of constraint algorithm
+ constraint-algorithm     = Lincs
+ ; Do not constrain the start configuration
+ continuation             = no
+ ; Use successive overrelaxation to reduce the number of shake iterations
+ Shake-SOR                = no
+ ; Relative tolerance of shake
+ shake-tol                = 0.0001
+ ; Highest order in the expansion of the constraint coupling matrix
+ lincs-order              = 4
+ ; Number of iterations in the final step of LINCS. 1 is fine for
+ ; normal simulations, but use 2 to conserve energy in NVE runs.
+ ; For energy minimization with constraints it should be 4 to 8.
+ lincs-iter               = 1
+ ; Lincs will write a warning to the stderr if in one step a bond
+ ; rotates over more degrees than
+ lincs-warnangle          = 30
+ ; Convert harmonic bonds to morse potentials
+ morse                    = no
+ ; ENERGY GROUP EXCLUSIONS
+ ; Pairs of energy groups for which all non-bonded interactions are excluded
+ energygrp-excl           = 
+ ; WALLS                
+ ; Number of walls, type, atom types, densities and box-z scale factor for Ewald
+ nwall                    = 0
+ wall-type                = 9-3
+ wall-r-linpot            = -1
+ wall-atomtype            = 
+ wall-density             = 
+ wall-ewald-zfac          = 3
+ ; COM PULLING          
+ pull                     = no
+ ; AWH biasing          
+ awh                      = no
+ ; ENFORCED ROTATION    
+ ; Enforced rotation: No or Yes
+ rotation                 = no
+ ; Group to display and/or manipulate in interactive MD session
+ IMD-group                = 
+ ; NMR refinement stuff 
+ ; Distance restraints type: No, Simple or Ensemble
+ disre                    = No
+ ; Force weighting of pairs in one distance restraint: Conservative or Equal
+ disre-weighting          = Conservative
+ ; Use sqrt of the time averaged times the instantaneous violation
+ disre-mixed              = no
+ disre-fc                 = 1000
+ disre-tau                = 0
+ ; Output frequency for pair distances to energy file
+ nstdisreout              = 100
+ ; Orientation restraints: No or Yes
+ orire                    = no
+ ; Orientation restraints force constant and tau for time averaging
+ orire-fc                 = 0
+ orire-tau                = 0
+ orire-fitgrp             = 
+ ; Output frequency for trace(SD) and S to energy file
+ nstorireout              = 100
+ ; Free energy variables
+ free-energy              = no
+ couple-moltype           = 
+ couple-lambda0           = vdw-q
+ couple-lambda1           = vdw-q
+ couple-intramol          = no
+ init-lambda              = -1
+ init-lambda-state        = -1
+ delta-lambda             = 0
+ nstdhdl                  = 50
+ fep-lambdas              = 
+ mass-lambdas             = 
+ coul-lambdas             = 
+ vdw-lambdas              = 
+ bonded-lambdas           = 
+ restraint-lambdas        = 
+ temperature-lambdas      = 
+ calc-lambda-neighbors    = 1
+ init-lambda-weights      = 
+ dhdl-print-energy        = no
+ sc-alpha                 = 0
+ sc-power                 = 1
+ sc-r-power               = 6
+ sc-sigma                 = 0.3
+ sc-coul                  = no
+ separate-dhdl-file       = yes
+ dhdl-derivatives         = yes
+ dh_hist_size             = 0
+ dh_hist_spacing          = 0.1
+ ; Non-equilibrium MD stuff
+ acc-grps                 = 
+ accelerate               = 
+ freezegrps               = 
+ freezedim                = 
+ cos-acceleration         = 0
+ deform                   = 
+ ; simulated tempering variables
+ simulated-tempering      = no
+ simulated-tempering-scaling = geometric
+ sim-temp-low             = 300
+ sim-temp-high            = 300
+ ; Ion/water position swapping for computational electrophysiology setups
+ ; Swap positions along direction: no, X, Y, Z
+ swapcoords               = no
+ adress                   = no
+ ; User defined thingies
+ user1-grps               = 
+ user2-grps               = 
+ userint1                 = 0
+ userint2                 = 0
+ userint3                 = 0
+ userint4                 = 0
+ userreal1                = 0
+ userreal2                = 0
+ userreal3                = 0
+ userreal4                = 0
+ ; Electric fields
+ ; Format for electric-field-x, etc. is: four real variables:
+ ; amplitude (V/nm), frequency omega (1/ps), time for the pulse peak (ps),
+ ; and sigma (ps) width of the pulse. Omega = 0 means static field,
+ ; sigma = 0 means no pulse, leaving the field to be a cosine function.
+ electric-field-x         = 0 0 0 0
+ electric-field-y         = 0 0 0 0
+ electric-field-z         = 0 0 0 0
+ </String>
+ </ReferenceData>
index 1cab82b2695e1726920a1e8ebdcbe6787d0eca3c,deff9cccd30f2696a485db2031503716f98cc53b..60896c1045bf14c35777ed6ccd978290bc47ae47
@@@ -549,23 -549,34 +549,34 @@@ static void do_update_md(in
  
      if (doNoseHoover || doPROffDiagonal || doAcceleration)
      {
+         matrix stepM;
+         if (!doParrinelloRahman)
+         {
+             /* We should not apply PR scaling at this step */
+             clear_mat(stepM);
+         }
+         else
+         {
+             copy_mat(M, stepM);
+         }
          if (!doAcceleration)
          {
              updateMDLeapfrogGeneral<AccelerationType::none>
                  (start, nrend, doNoseHoover, dt, dtPressureCouple,
-                 ir, md, ekind, box, x, xprime, v, f, nh_vxi, M);
+                 ir, md, ekind, box, x, xprime, v, f, nh_vxi, stepM);
          }
          else if (ekind->bNEMD)
          {
              updateMDLeapfrogGeneral<AccelerationType::group>
                  (start, nrend, doNoseHoover, dt, dtPressureCouple,
-                 ir, md, ekind, box, x, xprime, v, f, nh_vxi, M);
+                 ir, md, ekind, box, x, xprime, v, f, nh_vxi, stepM);
          }
          else
          {
              updateMDLeapfrogGeneral<AccelerationType::cosine>
                  (start, nrend, doNoseHoover, dt, dtPressureCouple,
-                 ir, md, ekind, box, x, xprime, v, f, nh_vxi, M);
+                 ir, md, ekind, box, x, xprime, v, f, nh_vxi, stepM);
          }
      }
      else
@@@ -1553,7 -1564,6 +1564,7 @@@ void update_constraints(FIL
                          t_idef                        *idef,
                          tensor                         vir_part,
                          t_commrec                     *cr,
 +                        const gmx_multisim_t          *ms,
                          t_nrnb                        *nrnb,
                          gmx_wallcycle_t                wcycle,
                          gmx_update_t                  *upd,
          if (EI_VV(inputrec->eI) && bFirstHalf)
          {
              constrain(nullptr, bLog, bEner, constr, idef,
 -                      inputrec, cr, step, 1, 1.0, md,
 +                      inputrec, cr, ms, step, 1, 1.0, md,
                        as_rvec_array(state->x.data()), as_rvec_array(state->v.data()), as_rvec_array(state->v.data()),
                        bMolPBC, state->box,
                        state->lambda[efptBONDED], dvdlambda,
          else
          {
              constrain(nullptr, bLog, bEner, constr, idef,
 -                      inputrec, cr, step, 1, 1.0, md,
 +                      inputrec, cr, ms, step, 1, 1.0, md,
                        as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()), nullptr,
                        bMolPBC, state->box,
                        state->lambda[efptBONDED], dvdlambda,
              wallcycle_start(wcycle, ewcCONSTR);
  
              constrain(nullptr, bLog, bEner, constr, idef,
 -                      inputrec, cr, step, 1, 0.5, md,
 +                      inputrec, cr, ms, step, 1, 0.5, md,
                        as_rvec_array(state->x.data()), as_rvec_array(upd->xp.data()), nullptr,
                        bMolPBC, state->box,
                        state->lambda[efptBONDED], dvdlambda,
index c8d8b54c4c9eda1825f182f441e61f364640e10c,9c6f463ae595f159d0ca5175fc26529796972cca..67824e3258e4705a9d4c2312dab72175a0124469
@@@ -140,17 -140,6 +140,17 @@@ tMPI_Thread_mutex_t deform_init_box_mut
  namespace gmx
  {
  
 +/*! \brief Barrier for safe simultaneous thread access to mdrunner data
 + *
 + * Used to ensure that the master thread does not modify mdrunner during copy
 + * on the spawned threads. */
 +static void threadMpiMdrunnerAccessBarrier()
 +{
 +#if GMX_THREAD_MPI
 +    MPI_Barrier(MPI_COMM_WORLD);
 +#endif
 +}
 +
  void Mdrunner::reinitializeOnSpawnedThread()
  {
      // TODO This duplication is formally necessary if any thread might
      // Mdrunner.
      fnm = dup_tfn(nfile, fnm);
  
 -    cr  = reinitialize_commrec_for_this_thread(cr);
 +    threadMpiMdrunnerAccessBarrier();
  
 -    if (!MASTER(cr))
 -    {
 -        // Only the master rank writes to the log files
 -        fplog = nullptr;
 -    }
 +    cr  = reinitialize_commrec_for_this_thread(cr, ms);
 +
 +    GMX_RELEASE_ASSERT(!MASTER(cr), "reinitializeOnSpawnedThread should only be called on spawned threads");
 +
 +    // Only the master rank writes to the log file
 +    fplog = nullptr;
  }
  
  /*! \brief The callback used for running on spawned threads.
   * argument permitted to the thread-launch API call, copies it to make
   * a new runner for this thread, reinitializes necessary data, and
   * proceeds to the simulation. */
 -static void mdrunner_start_fn(void *arg)
 +static void mdrunner_start_fn(const void *arg)
  {
      try
      {
          auto masterMdrunner = reinterpret_cast<const gmx::Mdrunner *>(arg);
          /* copy the arg list to make sure that it's thread-local. This
 -           doesn't copy pointed-to items, of course, but those are all
 -           const. */
 +           doesn't copy pointed-to items, of course; fnm, cr and fplog
 +           are reset in the call below, all others should be const. */
          gmx::Mdrunner mdrunner = *masterMdrunner;
          mdrunner.reinitializeOnSpawnedThread();
          mdrunner.mdrunner();
   * (including the main thread) for thread-parallel runs. This in turn
   * calls mdrunner() for each thread. All options are the same as for
   * mdrunner(). */
 -t_commrec *Mdrunner::spawnThreads(int numThreadsToLaunch)
 +t_commrec *Mdrunner::spawnThreads(int numThreadsToLaunch) const
  {
  
      /* first check whether we even need to start tMPI */
          return cr;
      }
  
 -    gmx::Mdrunner spawnedMdrunner = *this;
 -    // TODO This duplication is formally necessary if any thread might
 -    // modify any memory in fnm or the pointers it contains. If the
 -    // contents are ever provably const, then we can remove this
 -    // allocation (and memory leak).
 -    // TODO This should probably become part of a copy constructor for
 -    // Mdrunner.
 -    spawnedMdrunner.fnm = dup_tfn(this->nfile, fnm);
 -
  #if GMX_THREAD_MPI
      /* now spawn new threads that start mdrunner_start_fn(), while
         the main thread returns, we set thread affinity later */
      if (tMPI_Init_fn(TRUE, numThreadsToLaunch, TMPI_AFFINITY_NONE,
 -                     mdrunner_start_fn, static_cast<void*>(&spawnedMdrunner)) != TMPI_SUCCESS)
 +                     mdrunner_start_fn, static_cast<const void*>(this)) != TMPI_SUCCESS)
      {
          GMX_THROW(gmx::InternalError("Failed to spawn thread-MPI threads"));
      }
 +
 +    threadMpiMdrunnerAccessBarrier();
  #else
      GMX_UNUSED_VALUE(mdrunner_start_fn);
  #endif
  
 -    return reinitialize_commrec_for_this_thread(cr);
 +    return reinitialize_commrec_for_this_thread(cr, ms);
  }
  
  }      // namespace
@@@ -449,9 -444,9 +449,9 @@@ int Mdrunner::mdrunner(
  
      /* CAUTION: threads may be started later on in this function, so
         cr doesn't reflect the final parallel state right now */
 -    gmx::MDModules mdModules;
 -    t_inputrec     inputrecInstance;
 -    t_inputrec    *inputrec = &inputrecInstance;
 +    std::unique_ptr<gmx::MDModules> mdModules(new gmx::MDModules);
 +    t_inputrec                      inputrecInstance;
 +    t_inputrec                     *inputrec = &inputrecInstance;
      snew(mtop, 1);
  
      if (mdrunOptions.continuationOptions.appendFiles)
      gmx::LoggerOwner logOwner(buildLogger(fplog, cr));
      gmx::MDLogger    mdlog(logOwner.logger());
  
 -    hwinfo = gmx_detect_hardware(mdlog, cr);
 +    hwinfo = gmx_detect_hardware(mdlog);
  
 -    gmx_print_detected_hardware(fplog, cr, mdlog, hwinfo);
 +    gmx_print_detected_hardware(fplog, cr, ms, mdlog, hwinfo);
  
      std::vector<int> gpuIdsToUse;
      auto             compatibleGpus = getCompatibleGpus(hwinfo->gpu_info);
      GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
  
      // TODO: Error handling
 -    mdModules.assignOptionsToModules(*inputrec->params, nullptr);
 +    mdModules->assignOptionsToModules(*inputrec->params, nullptr);
  
      if (fplog != nullptr)
      {
      snew(fcd, 1);
  
      /* This needs to be called before read_checkpoint to extend the state */
 -    init_disres(fplog, mtop, inputrec, cr, fcd, globalState.get(), replExParams.exchangeInterval > 0);
 +    init_disres(fplog, mtop, inputrec, cr, ms, fcd, globalState.get(), replExParams.exchangeInterval > 0);
  
 -    init_orires(fplog, mtop, inputrec, cr, globalState.get(), &(fcd->orires));
 +    init_orires(fplog, mtop, inputrec, cr, ms, globalState.get(), &(fcd->orires));
  
      if (inputrecDeform(inputrec))
      {
  
      /* Initialize per-physical-node MPI process/thread ID and counters. */
      gmx_init_intranode_counters(cr);
 -    if (cr->ms && cr->ms->nsim > 1 && !opt2bSet("-multidir", nfile, fnm))
 -    {
 -        GMX_LOG(mdlog.info).asParagraph().
 -            appendText("The -multi flag is deprecated, and may be removed in a future version. Please "
 -                       "update your workflows to use -multidir instead.");
 -    }
  #if GMX_MPI
 -    if (MULTISIM(cr))
 +    if (isMultiSim(ms))
      {
          GMX_LOG(mdlog.warning).asParagraph().appendTextFormatted(
                  "This is simulation %d out of %d running as a composite GROMACS\n"
                  "multi-simulation job. Setup for this simulation:\n",
 -                cr->ms->sim, cr->ms->nsim);
 +                ms->sim, ms->nsim);
      }
      GMX_LOG(mdlog.warning).appendTextFormatted(
              "Using %d MPI %s\n",
      check_and_update_hw_opt_2(&hw_opt, inputrec->cutoff_scheme);
  
      /* Check and update the number of OpenMP threads requested */
 -    checkAndUpdateRequestedNumOpenmpThreads(&hw_opt, *hwinfo, cr, pmeRunMode, *mtop);
 +    checkAndUpdateRequestedNumOpenmpThreads(&hw_opt, *hwinfo, cr, ms, pmeRunMode, *mtop);
  
      gmx_omp_nthreads_init(mdlog, cr,
                            hwinfo->nthreads_hw_avail,
      {
          // Produce the task assignment for this rank.
          gpuTaskAssignment = runTaskAssignment(gpuIdsToUse, userGpuTaskAssignment, *hwinfo,
 -                                              mdlog, cr, gpuTasksOnThisRank);
 +                                              mdlog, cr, ms, gpuTasksOnThisRank);
      }
      GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
  
      {
          MPI_Barrier(cr->mpi_comm_mysim);
      }
 -    if (MULTISIM(cr))
 +    if (isMultiSim(ms))
      {
-         MPI_Barrier(ms->mpi_comm_masters);
+         if (MASTER(cr))
+         {
 -            MPI_Barrier(cr->ms->mpi_comm_masters);
++            MPI_Barrier(ms->mpi_comm_masters);
+         }
+         /* We need another barrier to prevent non-master ranks from contiuing
+          * when an error occured in a different simulation.
+          */
+         MPI_Barrier(cr->mpi_comm_mysim);
      }
  #endif
  
  
      checkHardwareOversubscription(numThreadsOnThisRank,
                                    *hwinfo->hardwareTopology,
 -                                  cr, mdlog);
 +                                  cr, ms, mdlog);
  
      if (hw_opt.thread_affinity != threadaffOFF)
      {
          gmx_check_thread_affinity_set(mdlog, cr,
                                        &hw_opt, hwinfo->nthreads_hw_avail, TRUE);
  
 +        int numThreadsOnThisNode, intraNodeThreadOffset;
 +        analyzeThreadsOnThisNode(cr, ms, nullptr, numThreadsOnThisRank, &numThreadsOnThisNode,
 +                                 &intraNodeThreadOffset);
 +
          /* Set the CPU affinity */
          gmx_set_thread_affinity(mdlog, cr, &hw_opt, *hwinfo->hardwareTopology,
 -                                numThreadsOnThisRank, nullptr);
 +                                numThreadsOnThisRank, numThreadsOnThisNode,
 +                                intraNodeThreadOffset, nullptr);
      }
  
      if (mdrunOptions.timingOptions.resetStep > -1)
      {
          /* Initiate forcerecord */
          fr                 = mk_forcerec();
 -        fr->forceProviders = mdModules.initForceProviders();
 +        fr->forceProviders = mdModules->initForceProviders();
          init_forcerec(fplog, mdlog, fr, fcd,
                        inputrec, mtop, cr, box,
                        opt2fn("-table", nfile, fnm),
          }
  
          /* Now do whatever the user wants us to do (how flexible...) */
 -        my_integrator(inputrec->eI) (fplog, cr, mdlog, nfile, fnm,
 +        my_integrator(inputrec->eI) (fplog, cr, ms, mdlog, nfile, fnm,
                                       oenv,
                                       mdrunOptions,
                                       vsite, constr,
 -                                     mdModules.outputProvider(),
 +                                     mdModules->outputProvider(),
                                       inputrec, mtop,
                                       fcd,
                                       globalState.get(),
                 inputrec, nrnb, wcycle, walltime_accounting,
                 fr ? fr->nbv : nullptr,
                 pmedata,
 -               EI_DYNAMICS(inputrec->eI) && !MULTISIM(cr));
 +               EI_DYNAMICS(inputrec->eI) && !isMultiSim(ms));
  
      // Free PME data
      if (pmedata)
      // As soon as we destroy GPU contexts after mdrunner() exits, these lines should go.
      mdAtoms.reset(nullptr);
      globalState.reset(nullptr);
 +    mdModules.reset(nullptr);   // destruct force providers here as they might also use the GPU
  
      /* Free GPU memory and set a physical node tMPI barrier (which should eventually go away) */
 -    free_gpu_resources(fr, cr);
 +    free_gpu_resources(fr, cr, ms);
      free_gpu(nonbondedDeviceInfo);
      free_gpu(pmeDeviceInfo);
  
      if (MASTER(cr) && continuationOptions.appendFiles)
      {
          gmx_log_close(fplog);
 +        fplog = nullptr;
      }
  
      rc = (int)gmx_get_stop_condition();
         wait for that. */
      if (PAR(cr) && MASTER(cr))
      {
 +        done_commrec(cr);
          tMPI_Finalize();
      }
  #endif