Merge branch release-2018
authorMark Abraham <mark.j.abraham@gmail.com>
Wed, 28 Mar 2018 11:19:26 +0000 (13:19 +0200)
committerMark Abraham <mark.j.abraham@gmail.com>
Wed, 28 Mar 2018 11:28:03 +0000 (13:28 +0200)
Some meaningful conflicts in init_multisim associated with multisim
fixes and clarifications in release-2018 combining with reorganization
in master branch.

Preserved the master-branch pre-submit matrix, since all branches have
separate changes that cope with gpu_id=2 no longer being valid.

Change-Id: Idefe4f1ec3bbbc76dfedf03a81abef8233328de2

1  2 
cmake/gmxVersionInfo.cmake
docs/CMakeLists.txt
docs/user-guide/mdrun-performance.rst
src/gromacs/fileio/checkpoint.cpp
src/gromacs/mdlib/main.cpp
src/gromacs/pulling/pull.cpp
src/programs/mdrun/md.cpp

index 56c9b90712c71cf5c5b6d909298fa2f59aff331f,d87c02da7ff5c776036c32b4934b70475774a4e4..109ebcbf05e224c73deb3921cd95feadfa1da709
@@@ -56,7 -56,6 +56,7 @@@
  #         GROMACS     5.1    1
  #         GROMACS     2016   2
  #         GROMACS     2018   3
 +#         GROMACS     2019   4
  #   LIBRARY_SOVERSION_MINOR so minor version for the built libraries.
  #       Should be increased for each release that changes only the implementation.
  #       In GROMACS, the typical policy is to increase it for each patch version
  # header) can provide useful information for, e.g., diagnosing bug reports and
  # identifying what exact version the user was using.  The following formats are
  # possible (with examples given for a particular version):
 -#   4.6.1       Plain version number without any suffix signifies a build from
 -#               a released source tarball.
 -#   4.6.1-dev   '-dev' suffix signifies all other builds. If there is no other
 -#               information, either the user built the code outside any git
 -#               repository, or disabled the version info generation.
 -#   4.6.1-dev-YYYYMMDD-1234abc
 -#               The YYYYMMDD part shows the commit date (not author date) of
 -#               the HEAD commit from which the code was built.  The abbreviated
 -#               hash is the hash of that commit (the full hash is available in
 -#               'gmx -version' output).
 -#               If the HEAD hash is not identified as coming from branches in
 -#               "authoritative" GROMACS repositories, 'gmx -version' will show
 -#               the nearest ancestor commit that is identified as such (but see
 -#               the '-local' and '-unknown' suffixes below).
 -#   4.6.1-dev-YYYYMMDD-1234abc-dirty
 -#               As above, but there were local modifications in the source tree
 -#               when the code was built.
 -#   4.6.1-dev-YYYYMMDD-1234abc-unknown
 -#               As above, but there were no remotes in the repository that
 -#               could be identified as "authoritative" GROMACS repositories.
 -#               This happens if the code is not cloned from git.gromacs.org
 -#               or gerrit.gromacs.org.
 -#   4.6.1-dev-YYYYMMDD-1234abc-local
 -#               As above, but there were no commits in the recent history of
 -#               the branch that could be identified as coming from
 -#               "authoritative" GROMACS repositories.  This should be
 -#               relatively rare.
 +#   2018.1       Plain version number without any suffix signifies a build from
 +#                a released source tarball.
 +#   2018.1-dev   '-dev' suffix signifies all other builds. If there is no other
 +#                information, either the user built the code outside any git
 +#                repository, or disabled the version info generation.
 +#   2018.1-dev-YYYYMMDD-1234abc
 +#                The YYYYMMDD part shows the commit date (not author date) of
 +#                the HEAD commit from which the code was built.  The abbreviated
 +#                hash is the hash of that commit (the full hash is available in
 +#                'gmx -version' output).
 +#                If the HEAD hash is not identified as coming from branches in
 +#                "authoritative" GROMACS repositories, 'gmx -version' will show
 +#                the nearest ancestor commit that is identified as such (but see
 +#                the '-local' and '-unknown' suffixes below).
 +#   2018.1-dev-YYYYMMDD-1234abc-dirty
 +#                As above, but there were local modifications in the source tree
 +#                when the code was built.
 +#   2018.1-dev-YYYYMMDD-1234abc-unknown
 +#                As above, but there were no remotes in the repository that
 +#                could be identified as "authoritative" GROMACS repositories.
 +#                This happens if the code is not cloned from git.gromacs.org
 +#                or gerrit.gromacs.org.
 +#   2018.1-dev-YYYYMMDD-1234abc-local
 +#                As above, but there were no commits in the recent history of
 +#                the branch that could be identified as coming from
 +#                "authoritative" GROMACS repositories.  This should be
 +#                relatively rare.
  #
  # Other variables set here are not intended for use outside this file.
  # The scripts gmxGenerateVersionInfo.cmake and gmxConfigureVersionInfo.cmake
  
  # The GROMACS convention is that these are the version number of the next
  # release that is going to be made from this branch.
 -set(GMX_VERSION_MAJOR 2018)
 -set(GMX_VERSION_PATCH 2)
 +set(GMX_VERSION_MAJOR 2019)
 +set(GMX_VERSION_PATCH 0)
  # The suffix, on the other hand, is used mainly for betas and release
  # candidates, where it signifies the most recent such release from
  # this branch; it will be empty before the first such release, as well
@@@ -203,8 -202,8 +203,8 @@@ set(GMX_VERSION_SUFFIX ""
  # here. The important thing is to minimize the chance of third-party
  # code being able to dynamically link with a version of libgromacs
  # that might not work.
 -set(LIBRARY_SOVERSION_MAJOR 3)
 -set(LIBRARY_SOVERSION_MINOR 1)
 +set(LIBRARY_SOVERSION_MAJOR 4)
 +set(LIBRARY_SOVERSION_MINOR 0)
  set(LIBRARY_VERSION ${LIBRARY_SOVERSION_MAJOR}.${LIBRARY_SOVERSION_MINOR}.0)
  
  #####################################################################
@@@ -227,7 -226,7 +227,7 @@@ endif(
  
  set(REGRESSIONTEST_VERSION "${GMX_VERSION_STRING}")
  set(REGRESSIONTEST_BRANCH "refs/heads/release-2018")
- set(REGRESSIONTEST_MD5SUM "5a609bea6b6777072063be9e665121ef" CACHE INTERNAL "MD5 sum of the regressiontests tarball")
+ set(REGRESSIONTEST_MD5SUM "1a94916e2cf90e34fddb3514a65e0154" CACHE INTERNAL "MD5 sum of the regressiontests tarball")
  
  math(EXPR GMX_VERSION_NUMERIC
       "${GMX_VERSION_MAJOR}*10000 + ${GMX_VERSION_PATCH}")
diff --combined docs/CMakeLists.txt
index f18d56fd4a1cab6158aaef5a8ccb555d9b7e7135,3909260f721db62a10001db78af7cb47ce570d86..845ca619ddef1d42b84636ba56f019976d534836
@@@ -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
          fragments/doxygen-links.rst
          install-guide/index.rst
          release-notes/index.rst
+         release-notes/2018/2018.2.rst
          release-notes/2018/2018.1.rst
          release-notes/2018/major/highlights.rst
          release-notes/2018/major/features.rst
index 859fc5f0680d338442ee9ddaf5d577912ea3d383,9d09a1199b59632d6a3ce46bb2002103dc20bb39..3fa3c280b06643ea9b8208321c2a0caf0b2fe384
@@@ -1,7 -1,7 +1,7 @@@
  .. _gmx-performance:
  
- Getting good performance from mdrun
- ===================================
+ Getting good performance from :ref:`mdrun <gmx mdrun>`
+ ======================================================
  The |Gromacs| build system and the :ref:`gmx mdrun` tool has a lot of built-in
  and configurable intelligence to detect your hardware and make pretty
  effective use of that hardware. For a lot of casual and serious use of
@@@ -137,7 -137,7 +137,7 @@@ see the Reference Manual. The most impo
          members of its domain. A GPU may perform work for more than
          one PP rank, but it is normally most efficient to use a single
          PP rank per GPU and for that rank to have thousands of
-         particles. When the work of a PP rank is done on the CPU, mdrun
+         particles. When the work of a PP rank is done on the CPU, :ref:`mdrun <gmx mdrun>`
          will make extensive use of the SIMD capabilities of the
          core. There are various `command-line options
          <controlling-the-domain-decomposition-algorithm` to control
          there are separate PME ranks, then the remaining ranks handle
          the PP work. Otherwise, all ranks do both PP and PME work.
  
- Running mdrun within a single node
- ----------------------------------
+ Running :ref:`mdrun <gmx mdrun>` within a single node
+ -----------------------------------------------------
  
  :ref:`gmx mdrun` can be configured and compiled in several different ways that
  are efficient to use within a single :term:`node`. The default configuration
@@@ -193,7 -193,7 +193,7 @@@ behavior
  ``-ntomp``
      The total number of OpenMP threads per rank to start. The
      default, 0, will start one thread on each available core.
-     Alternatively, mdrun will honor the appropriate system
+     Alternatively, :ref:`mdrun <gmx mdrun>` will honor the appropriate system
      environment variable (e.g. ``OMP_NUM_THREADS``) if set.
  
  ``-npme``
  
  ``-pin``
      Can be set to "auto," "on" or "off" to control whether
-     mdrun will attempt to set the affinity of threads to cores.
-     Defaults to "auto," which means that if mdrun detects that all the
-     cores on the node are being used for mdrun, then it should behave
+     :ref:`mdrun <gmx mdrun>` will attempt to set the affinity of threads to cores.
+     Defaults to "auto," which means that if :ref:`mdrun <gmx mdrun>` detects that all the
+     cores on the node are being used for :ref:`mdrun <gmx mdrun>`, then it should behave
      like "on," and attempt to set the affinities (unless they are
      already set by something else).
  
  ``-pinoffset``
      If ``-pin on``, specifies the logical core number to
-     which mdrun should pin the first thread. When running more than
-     one instance of mdrun on a node, use this option to to avoid
-     pinning threads from different mdrun instances to the same core.
+     which :ref:`mdrun <gmx mdrun>` should pin the first thread. When running more than
+     one instance of :ref:`mdrun <gmx mdrun>` on a node, use this option to to avoid
+     pinning threads from different :ref:`mdrun <gmx mdrun>` instances to the same core.
  
  ``-pinstride``
      If ``-pin on``, specifies the stride in logical core
-     numbers for the cores to which mdrun should pin its threads. When
-     running more than one instance of mdrun on a node, use this option
-     to to avoid pinning threads from different mdrun instances to the
+     numbers for the cores to which :ref:`mdrun <gmx mdrun>` should pin its threads. When
+     running more than one instance of :ref:`mdrun <gmx mdrun>` on a node, use this option
+     to to avoid pinning threads from different :ref:`mdrun <gmx mdrun>` instances to the
      same core.  Use the default, 0, to minimize the number of threads
-     per physical core - this lets mdrun manage the hardware-, OS- and
+     per physical core - this lets :ref:`mdrun <gmx mdrun>` manage the hardware-, OS- and
      configuration-specific details of how to map logical cores to
      physical cores.
  
      A string that specifies the ID numbers of the GPUs that
      are available to be used by ranks on this node. For example,
      "12" specifies that the GPUs with IDs 1 and 2 (as reported
-     by the GPU runtime) can be used by mdrun. This is useful
+     by the GPU runtime) can be used by :ref:`mdrun <gmx mdrun>`. This is useful
      when sharing a node with other computations, or if a GPU
-     is best used to support a display. If many GPUs are
+     is best used to support a display.  Without specifying this
+     parameter, :ref:`mdrun <gmx mdrun>` will utilize all GPUs. When many GPUs are
      present, a comma may be used to separate the IDs, so
-     "12,13" would make GPUs 12 and 13 available to mdrun.
+     "12,13" would make GPUs 12 and 13 available to :ref:`mdrun <gmx mdrun>`.
      It could be necessary to use different GPUs on different
      nodes of a simulation, in which case the environment
      variable ``GMX_GPU_ID`` can be set differently for the ranks
      on different nodes to achieve that result.
+     In |Gromacs| versions preceding 2018 this parameter used to
+     specify both GPU availability and GPU task assignment.
+     The latter is now done with the ``-gputasks`` parameter.
  
  ``-gputasks``
      A string that specifies the ID numbers of the GPUs to be
      used by corresponding GPU tasks on this node. For example,
      "0011" specifies that the first two GPU tasks will use GPU 0,
      and the other two use GPU 1. When using this option, the
-     number of ranks must be known to mdrun, as well as where
+     number of ranks must be known to :ref:`mdrun <gmx mdrun>`, as well as where
      tasks of different types should be run, such as by using
-     ``-nb gpu``.
+     ``-nb gpu`` - only the tasks which are set to run on GPUs
+     count for parsing the mapping.
+     In |Gromacs| versions preceding 2018 only a single type
+     of GPU task could be run on any rank. Now that there is some
+     support for running PME on GPUs, the number of GPU tasks
+     (and the number of GPU IDs expected in the ``-gputasks`` string)
+     can actually be 2 for a single-rank simulation. The IDs
+     still have to be the same in this case, as using multiple GPUs
+     per single rank is not yet implemented.
+     The order of GPU tasks per rank in the string is short-range first,
+     PME second. The order of ranks with different kinds of GPU tasks
+     is the same by default, but can be influenced with the ``-ddorder``
+     option and gets quite complex when using multiple nodes.
+     The GPU task assignment (whether manually set, or automated),
+     will be reported in the :ref:`mdrun <gmx mdrun>` output on
+     the first physical node of the simulation. For example:
  
- Examples for mdrun on one node
- ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+     ::
+       gmx mdrun -gputasks 0001 -nb gpu -pme gpu -npme 1 -ntmpi 4
+     will produce the following output in the log file/terminal:
+     ::
+       On host tcbl14 2 GPUs user-selected for this run.
+       Mapping of GPU IDs to the 4 GPU tasks in the 4 ranks on this node:
+       PP:0,PP:0,PP:0,PME:1
+     In this case, 3 ranks are set by user to compute short-range work
+     on GPU 0, and 1 rank to compute PME on GPU 1.
+     The detailed indexing of the GPUs is also reported in the log file.
+     For more information about GPU tasks, please refer to
+     :ref:`Types of GPU tasks<gmx-gpu-tasks>`.
+ Examples for :ref:`mdrun <gmx mdrun>` on one node
+ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  
  ::
  
      gmx mdrun
  
- Starts mdrun using all the available resources. mdrun
+ Starts :ref:`mdrun <gmx mdrun>` using all the available resources. :ref:`mdrun <gmx mdrun>`
  will automatically choose a fairly efficient division
  into thread-MPI ranks, OpenMP threads and assign work
  to compatible GPUs. Details will vary with hardware
@@@ -289,7 -327,7 +327,7 @@@ and the kind of simulation being run
  
      gmx mdrun -nt 8
  
- Starts mdrun using 8 threads, which might be thread-MPI
+ Starts :ref:`mdrun <gmx mdrun>` using 8 threads, which might be thread-MPI
  or OpenMP threads depending on hardware and the kind
  of simulation being run.
  
  
      gmx mdrun -ntmpi 2 -ntomp 4
  
- Starts mdrun using eight total threads, with four thread-MPI
+ Starts :ref:`mdrun <gmx mdrun>` using eight total threads, with four thread-MPI
  ranks and two OpenMP threads per core. You should only use
  these options when seeking optimal performance, and
  must take care that the ranks you create can have
@@@ -310,7 -348,7 +348,7 @@@ a multiple of the number of threads pe
  
      gmx mdrun -gpu_id 12
  
- Starts mdrun using GPUs with IDs 1 and 2 (e.g. because
+ Starts :ref:`mdrun <gmx mdrun>` using GPUs with IDs 1 and 2 (e.g. because
  GPU 0 is dedicated to running a display). This requires
  two thread-MPI ranks, and will split the available
  CPU cores between them using OpenMP threads.
  
      gmx mdrun -ntmpi 4 -nb gpu -gputasks 1122
  
- Starts mdrun using four thread-MPI ranks, and maps them
+ Starts :ref:`mdrun <gmx mdrun>` using four thread-MPI ranks, and maps them
  to GPUs with IDs 1 and 2. The CPU cores available will
  be split evenly between the ranks using OpenMP threads.
  
      gmx mdrun -nt 6 -pin on -pinoffset 0
      gmx mdrun -nt 6 -pin on -pinoffset 3
  
- Starts two mdrun processes, each with six total threads.
+ Starts two :ref:`mdrun <gmx mdrun>` processes, each with six total threads.
  Threads will have their affinities set to particular
  logical cores, beginning from the logical core
  with rank 0 or 3, respectively. The above would work
  well on an Intel CPU with six physical cores and
  hyper-threading enabled. Use this kind of setup only
- if restricting mdrun to a subset of cores to share a
+ if restricting :ref:`mdrun <gmx mdrun>` to a subset of cores to share a
  node with other processes.
  
  ::
@@@ -347,12 -385,12 +385,12 @@@ as the hardware and MPI setup will perm
  MPI setup is restricted to one node, then the resulting
  :ref:`gmx mdrun` will be local to that node.
  
- Running mdrun on more than one node
- -----------------------------------
+ Running :ref:`mdrun <gmx mdrun>` on more than one node
+ ------------------------------------------------------
  This requires configuring |Gromacs| to build with an external MPI
- library. By default, this mdrun executable is run with
+ library. By default, this :ref:`mdrun <gmx mdrun>` executable is run with
  :ref:`mdrun_mpi`. All of the considerations for running single-node
mdrun still apply, except that ``-ntmpi`` and ``-nt`` cause a fatal
:ref:`mdrun <gmx mdrun>` still apply, except that ``-ntmpi`` and ``-nt`` cause a fatal
  error, and instead the number of ranks is controlled by the
  MPI environment.
  Settings such as ``-npme`` are much more important when
@@@ -373,7 -411,7 +411,7 @@@ cases
      Defaults to "on." If "on," a Verlet-scheme simulation will
      optimize various aspects of the PME and DD algorithms, shifting
      load between ranks and/or GPUs to maximize throughput. Some
-     mdrun features are not compatible with this, and these ignore
+     :ref:`mdrun <gmx mdrun>` features are not compatible with this, and these ignore
      this option.
  
  ``-dlb``
@@@ -419,9 -457,9 +457,9 @@@ It is only aware of the number of rank
  the MPI environment, and does not explicitly manage
  any aspect of OpenMP during the optimization.
  
- Examples for mdrun on more than one node
- ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
- The examples and explanations for for single-node mdrun are
+ Examples for :ref:`mdrun <gmx mdrun>` on more than one node
+ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+ The examples and explanations for for single-node :ref:`mdrun <gmx mdrun>` are
  still relevant, but ``-nt`` is no longer the way
  to choose the number of MPI ranks.
  
@@@ -491,7 -529,7 +529,7 @@@ across ranks each to one OpenMP thread
  suitable when there are ten nodes, each with two GPUs, but another
  job on each node is using GPU 0. The job scheduler should set the
  affinity of threads of both jobs to their allocated cores, or the
- performance of mdrun will suffer greatly.
+ performance of :ref:`mdrun <gmx mdrun>` will suffer greatly.
  
  ::
  
@@@ -551,8 -589,8 +589,8 @@@ parallel hardware
      of ``-dds`` might need to be adjusted to account for high or low
      spatial inhomogeneity of the system.
  
- Finding out how to run mdrun better
- -----------------------------------
+ Finding out how to run :ref:`mdrun <gmx mdrun>` better
+ ------------------------------------------------------
  
  The Wallcycle module is used for runtime performance measurement of :ref:`gmx mdrun`.
  At the end of the log file of each run, the "Real cycle and time accounting" section
@@@ -577,22 -615,29 +615,28 @@@ The performance counters are
  * Neighbor search
  * Launch GPU operations
  * Communication of coordinates
 -* Born radii
  * Force
  * Waiting + Communication of force
  * Particle mesh Ewald
  * PME redist. X/F
- * PME spread/gather
+ * PME spread
+ * PME gather
  * PME 3D-FFT
  * PME 3D-FFT Communication
  * PME solve Lennard-Jones
+ * PME solve LJ
  * PME solve Elec
  * PME wait for particle-particle
  * Wait + Receive PME force
  * Wait GPU nonlocal
  * Wait GPU local
+ * Wait PME GPU spread
+ * Wait PME GPU gather
+ * Reduce PME GPU Force
  * Non-bonded position/force buffer operations
  * Virtual site spread
  * COM pull force
+ * AWH (accelerated weight histogram method)
  * Write trajectory
  * Update
  * Constraints
@@@ -632,7 -677,10 +676,10 @@@ An additional set of subcounters can of
  * Bonded-FEP force
  * Restraints force
  * Listed buffer operations
+ * Nonbonded pruning
  * Nonbonded force
+ * Launch non-bonded GPU tasks
+ * Launch PME GPU tasks
  * Ewald force correction
  * Non-bonded position buffer operations
  * Non-bonded force buffer operations
@@@ -649,8 -697,8 +696,8 @@@ maybe elsewher
  
  .. _gmx-mdrun-on-gpu:
  
- Running mdrun with GPUs
- -----------------------
+ Running :ref:`mdrun <gmx mdrun>` with GPUs
+ ------------------------------------------
  
  NVIDIA GPUs from the professional line (Tesla or Quadro) starting with
  the Kepler generation (compute capability 3.5 and later) support changing the
@@@ -745,20 -793,21 +792,21 @@@ of the short range interactions on the 
  Known limitations
  .................
  
- **Please note again the limitations outlined above!**
+ **Please note again the limitations outlined below!**
  
  - Only compilation with CUDA is supported.
  
- - Only a PME order of 4 is supported in GPU.
+ - Only a PME order of 4 is supported on GPUs.
  
  - PME will run on a GPU only when exactly one rank has a
    PME task, ie. decompositions with multiple ranks doing PME are not supported.
  
  - Only single precision is supported.
  
- - Free energy calculations are not supported, because only single PME grids can be calculated.
+ - Free energy calculations where charges are perturbed are not supported,
+   because only single PME grids can be calculated.
  
- - LJ PME is not supported on GPU.
+ - LJ PME is not supported on GPUs.
  
  Assigning tasks to GPUs
  .......................
index 9db7e551445dad36feb1b8e19af4bfd42e5e1e26,6714f88961c0fb2f48e46a5ce3453406c763fbad..fc7c7c4c1aa9b80167ab8acbc400ff3fbd441b1d
@@@ -51,8 -51,6 +51,8 @@@
  #include <sys/locking.h>
  #endif
  
 +#include <array>
 +
  #include "buildinfo.h"
  #include "gromacs/fileio/filetypes.h"
  #include "gromacs/fileio/gmxfio.h"
@@@ -183,7 -181,7 +183,7 @@@ const char *eawhh_names[eawhhNR] 
      "awh_coordpoint", "awh_umbrellaGridpoint",
      "awh_updatelist",
      "awh_logScaledSampleWeight",
-     "awh_numupdates"
+     "awh_numupdates",
      "awh_forceCorrelationGrid"
  };
  
@@@ -231,16 -229,20 +231,16 @@@ static void cp_error(
      gmx_fatal(FARGS, "Checkpoint file corrupted/truncated, or maybe you are out of disk space?");
  }
  
 -static void do_cpt_string_err(XDR *xd, gmx_bool bRead, const char *desc, char **s, FILE *list)
 +static void do_cpt_string_err(XDR *xd, const char *desc, gmx::ArrayRef<char> s, FILE *list)
  {
 -    if (bRead)
 -    {
 -        snew(*s, CPTSTRLEN);
 -    }
 -    if (xdr_string(xd, s, CPTSTRLEN) == 0)
 +    char *data = s.data();
 +    if (xdr_string(xd, &data, s.size()) == 0)
      {
          cp_error();
      }
      if (list)
      {
 -        fprintf(list, "%s = %s\n", desc, *s);
 -        sfree(*s);
 +        fprintf(list, "%s = %s\n", desc, data);
      }
  }
  
@@@ -874,49 -876,23 +874,49 @@@ static int do_cpte_matrices(XDR *xd, St
      return ret;
  }
  
 -static void do_cpt_header(XDR *xd, gmx_bool bRead, int *file_version,
 -                          char **version, char **btime, char **buser, char **bhost,
 -                          int *double_prec,
 -                          char **fprog, char **ftime,
 -                          int *eIntegrator, int *simulation_part,
 -                          gmx_int64_t *step, double *t,
 -                          int *nnodes, int *dd_nc, int *npme,
 -                          int *natoms, int *ngtc, int *nnhpres, int *nhchainlength,
 -                          int *nlambda, int *flags_state,
 -                          int *flags_eks, int *flags_enh, int *flags_dfh, int *flags_awhh,
 -                          int *nED, int *eSwapCoords,
 -                          FILE *list)
 +// TODO Expand this into being a container of all data for
 +// serialization of a checkpoint, which can be stored by the caller
 +// (e.g. so that mdrun doesn't have to open the checkpoint twice).
 +// This will separate issues of allocation from those of
 +// serialization, help separate comparison from reading, and have
 +// better defined transformation functions to/from trajectory frame
 +// data structures.
 +struct CheckpointHeaderContents
 +{
 +    int         file_version;
 +    char        version[CPTSTRLEN];
 +    char        btime[CPTSTRLEN];
 +    char        buser[CPTSTRLEN];
 +    char        bhost[CPTSTRLEN];
 +    int         double_prec;
 +    char        fprog[CPTSTRLEN];
 +    char        ftime[CPTSTRLEN];
 +    int         eIntegrator;
 +    int         simulation_part;
 +    gmx_int64_t step;
 +    double      t;
 +    int         nnodes;
 +    ivec        dd_nc;
 +    int         npme;
 +    int         natoms;
 +    int         ngtc;
 +    int         nnhpres;
 +    int         nhchainlength;
 +    int         nlambda;
 +    int         flags_state;
 +    int         flags_eks;
 +    int         flags_enh;
 +    int         flags_dfh;
 +    int         flags_awhh;
 +    int         nED;
 +    int         eSwapCoords;
 +};
 +
 +static void do_cpt_header(XDR *xd, gmx_bool bRead, FILE *list,
 +                          CheckpointHeaderContents *contents)
  {
      bool_t res = 0;
      int    magic;
 -    int    idum = 0;
 -    char  *fhost;
  
      if (bRead)
      {
                    "The checkpoint file is corrupted or not a checkpoint file",
                    magic, CPT_MAGIC1);
      }
 +    char fhost[255];
      if (!bRead)
      {
 -        snew(fhost, 255);
          gmx_gethostname(fhost, 255);
      }
 -    do_cpt_string_err(xd, bRead, "GROMACS version", version, list);
 -    do_cpt_string_err(xd, bRead, "GROMACS build time", btime, list);
 -    do_cpt_string_err(xd, bRead, "GROMACS build user", buser, list);
 -    do_cpt_string_err(xd, bRead, "GROMACS build host", bhost, list);
 -    do_cpt_string_err(xd, bRead, "generating program", fprog, list);
 -    do_cpt_string_err(xd, bRead, "generation time", ftime, list);
 -    *file_version = cpt_version;
 -    do_cpt_int_err(xd, "checkpoint file version", file_version, list);
 -    if (*file_version > cpt_version)
 +    do_cpt_string_err(xd, "GROMACS version", contents->version, list);
 +    do_cpt_string_err(xd, "GROMACS build time", contents->btime, list);
 +    do_cpt_string_err(xd, "GROMACS build user", contents->buser, list);
 +    do_cpt_string_err(xd, "GROMACS build host", contents->bhost, list);
 +    do_cpt_string_err(xd, "generating program", contents->fprog, list);
 +    do_cpt_string_err(xd, "generation time", contents->ftime, list);
 +    contents->file_version = cpt_version;
 +    do_cpt_int_err(xd, "checkpoint file version", &contents->file_version, list);
 +    if (contents->file_version > cpt_version)
      {
 -        gmx_fatal(FARGS, "Attempting to read a checkpoint file of version %d with code of version %d\n", *file_version, cpt_version);
 +        gmx_fatal(FARGS, "Attempting to read a checkpoint file of version %d with code of version %d\n", contents->file_version, cpt_version);
      }
 -    if (*file_version >= 13)
 +    if (contents->file_version >= 13)
      {
 -        do_cpt_int_err(xd, "GROMACS double precision", double_prec, list);
 +        do_cpt_int_err(xd, "GROMACS double precision", &contents->double_prec, list);
      }
      else
      {
 -        *double_prec = -1;
 +        contents->double_prec = -1;
      }
 -    if (*file_version >= 12)
 +    if (contents->file_version >= 12)
      {
 -        do_cpt_string_err(xd, bRead, "generating host", &fhost, list);
 -        if (list == nullptr)
 -        {
 -            sfree(fhost);
 -        }
 +        do_cpt_string_err(xd, "generating host", fhost, list);
      }
 -    do_cpt_int_err(xd, "#atoms", natoms, list);
 -    do_cpt_int_err(xd, "#T-coupling groups", ngtc, list);
 -    if (*file_version >= 10)
 +    do_cpt_int_err(xd, "#atoms", &contents->natoms, list);
 +    do_cpt_int_err(xd, "#T-coupling groups", &contents->ngtc, list);
 +    if (contents->file_version >= 10)
      {
 -        do_cpt_int_err(xd, "#Nose-Hoover T-chains", nhchainlength, list);
 +        do_cpt_int_err(xd, "#Nose-Hoover T-chains", &contents->nhchainlength, list);
      }
      else
      {
 -        *nhchainlength = 1;
 +        contents->nhchainlength = 1;
      }
 -    if (*file_version >= 11)
 +    if (contents->file_version >= 11)
      {
 -        do_cpt_int_err(xd, "#Nose-Hoover T-chains for barostat ", nnhpres, list);
 +        do_cpt_int_err(xd, "#Nose-Hoover T-chains for barostat ", &contents->nnhpres, list);
      }
      else
      {
 -        *nnhpres = 0;
 +        contents->nnhpres = 0;
      }
 -    if (*file_version >= 14)
 +    if (contents->file_version >= 14)
      {
 -        do_cpt_int_err(xd, "# of total lambda states ", nlambda, list);
 +        do_cpt_int_err(xd, "# of total lambda states ", &contents->nlambda, list);
      }
      else
      {
 -        *nlambda = 0;
 +        contents->nlambda = 0;
      }
 -    do_cpt_int_err(xd, "integrator", eIntegrator, list);
 -    if (*file_version >= 3)
 +    do_cpt_int_err(xd, "integrator", &contents->eIntegrator, list);
 +    if (contents->file_version >= 3)
      {
 -        do_cpt_int_err(xd, "simulation part #", simulation_part, list);
 +        do_cpt_int_err(xd, "simulation part #", &contents->simulation_part, list);
      }
      else
      {
 -        *simulation_part = 1;
 +        contents->simulation_part = 1;
      }
 -    if (*file_version >= 5)
 +    if (contents->file_version >= 5)
      {
 -        do_cpt_step_err(xd, "step", step, list);
 +        do_cpt_step_err(xd, "step", &contents->step, list);
      }
      else
      {
 +        int idum = 0;
          do_cpt_int_err(xd, "step", &idum, list);
 -        *step = idum;
 +        contents->step = static_cast<gmx_int64_t>(idum);
      }
 -    do_cpt_double_err(xd, "t", t, list);
 -    do_cpt_int_err(xd, "#PP-ranks", nnodes, list);
 -    idum = 1;
 -    do_cpt_int_err(xd, "dd_nc[x]", dd_nc ? &(dd_nc[0]) : &idum, list);
 -    do_cpt_int_err(xd, "dd_nc[y]", dd_nc ? &(dd_nc[1]) : &idum, list);
 -    do_cpt_int_err(xd, "dd_nc[z]", dd_nc ? &(dd_nc[2]) : &idum, list);
 -    do_cpt_int_err(xd, "#PME-only ranks", npme, list);
 -    do_cpt_int_err(xd, "state flags", flags_state, list);
 -    if (*file_version >= 4)
 +    do_cpt_double_err(xd, "t", &contents->t, list);
 +    do_cpt_int_err(xd, "#PP-ranks", &contents->nnodes, list);
 +    do_cpt_int_err(xd, "dd_nc[x]", &contents->dd_nc[XX], list);
 +    do_cpt_int_err(xd, "dd_nc[y]", &contents->dd_nc[YY], list);
 +    do_cpt_int_err(xd, "dd_nc[z]", &contents->dd_nc[ZZ], list);
 +    do_cpt_int_err(xd, "#PME-only ranks", &contents->npme, list);
 +    do_cpt_int_err(xd, "state flags", &contents->flags_state, list);
 +    if (contents->file_version >= 4)
      {
 -        do_cpt_int_err(xd, "ekin data flags", flags_eks, list);
 -        do_cpt_int_err(xd, "energy history flags", flags_enh, list);
 +        do_cpt_int_err(xd, "ekin data flags", &contents->flags_eks, list);
 +        do_cpt_int_err(xd, "energy history flags", &contents->flags_enh, list);
      }
      else
      {
 -        *flags_eks   = 0;
 -        *flags_enh   = (*flags_state >> (estORIRE_DTAV+1));
 -        *flags_state = (*flags_state & ~((1<<(estORIRE_DTAV+1)) |
 -                                         (1<<(estORIRE_DTAV+2)) |
 -                                         (1<<(estORIRE_DTAV+3))));
 +        contents->flags_eks   = 0;
 +        contents->flags_enh   = (contents->flags_state >> (estORIRE_DTAV+1));
 +        contents->flags_state = (contents->flags_state & ~((1<<(estORIRE_DTAV+1)) |
 +                                                           (1<<(estORIRE_DTAV+2)) |
 +                                                           (1<<(estORIRE_DTAV+3))));
      }
 -    if (*file_version >= 14)
 +    if (contents->file_version >= 14)
      {
 -        do_cpt_int_err(xd, "df history flags", flags_dfh, list);
 +        do_cpt_int_err(xd, "df history flags", &contents->flags_dfh, list);
      }
      else
      {
 -        *flags_dfh = 0;
 +        contents->flags_dfh = 0;
      }
  
 -    if (*file_version >= 15)
 +    if (contents->file_version >= 15)
      {
 -        do_cpt_int_err(xd, "ED data sets", nED, list);
 +        do_cpt_int_err(xd, "ED data sets", &contents->nED, list);
      }
      else
      {
 -        *nED = 0;
 +        contents->nED = 0;
      }
  
 -    if (*file_version >= 16)
 +    if (contents->file_version >= 16)
      {
 -        do_cpt_int_err(xd, "swap", eSwapCoords, list);
 +        do_cpt_int_err(xd, "swap", &contents->eSwapCoords, list);
      }
      else
      {
 -        *eSwapCoords = eswapNO;
 +        contents->eSwapCoords = eswapNO;
      }
  
 -    if (*file_version >= 17)
 +    if (contents->file_version >= 17)
      {
 -        do_cpt_int_err(xd, "AWH history flags", flags_awhh, list);
 +        do_cpt_int_err(xd, "AWH history flags", &contents->flags_awhh, list);
      }
      else
      {
 -        *flags_awhh = 0;
 +        contents->flags_awhh = 0;
      }
  }
  
@@@ -1692,34 -1672,39 +1692,34 @@@ static int do_cpt_awh(XDR *xd, gmx_boo
  }
  
  static int do_cpt_files(XDR *xd, gmx_bool bRead,
 -                        gmx_file_position_t **p_outputfiles, int *nfiles,
 +                        std::vector<gmx_file_position_t> *outputfiles,
                          FILE *list, int file_version)
  {
 -    int                  i;
 -    gmx_off_t            offset;
 -    gmx_off_t            mask = 0xFFFFFFFFL;
 -    int                  offset_high, offset_low;
 -    char                *buf;
 -    gmx_file_position_t *outputfiles;
 +    gmx_off_t                   offset;
 +    gmx_off_t                   mask = 0xFFFFFFFFL;
 +    int                         offset_high, offset_low;
 +    std::array<char, CPTSTRLEN> buf;
 +    GMX_RELEASE_ASSERT(outputfiles, "Must have valid outputfiles");
  
 -    if (do_cpt_int(xd, "number of output files", nfiles, list) != 0)
 +    // Ensure that reading pre-allocates outputfiles, while writing
 +    // writes what is already there.
 +    int nfiles = outputfiles->size();
 +    if (do_cpt_int(xd, "number of output files", &nfiles, list) != 0)
      {
          return -1;
      }
 -
      if (bRead)
      {
 -        snew(*p_outputfiles, *nfiles);
 +        outputfiles->resize(nfiles);
      }
  
 -    outputfiles = *p_outputfiles;
 -
 -    for (i = 0; i < *nfiles; i++)
 +    for (auto &outputfile : *outputfiles)
      {
          /* 64-bit XDR numbers are not portable, so it is stored as separate high/low fractions */
          if (bRead)
          {
 -            do_cpt_string_err(xd, bRead, "output filename", &buf, list);
 -            std::strncpy(outputfiles[i].filename, buf, CPTSTRLEN-1);
 -            if (list == nullptr)
 -            {
 -                sfree(buf);
 -            }
 +            do_cpt_string_err(xd, "output filename", buf, list);
 +            std::strncpy(outputfile.filename, buf.data(), buf.size()-1);
  
              if (do_cpt_int(xd, "file_offset_high", &offset_high, list) != 0)
              {
              {
                  return -1;
              }
 -            outputfiles[i].offset = (static_cast<gmx_off_t>(offset_high) << 32 ) | ( static_cast<gmx_off_t>(offset_low) & mask );
 +            outputfile.offset = (static_cast<gmx_off_t>(offset_high) << 32 ) | ( static_cast<gmx_off_t>(offset_low) & mask );
          }
          else
          {
 -            buf = outputfiles[i].filename;
 -            do_cpt_string_err(xd, bRead, "output filename", &buf, list);
 +            do_cpt_string_err(xd, "output filename", outputfile.filename, list);
              /* writing */
 -            offset      = outputfiles[i].offset;
 +            offset      = outputfile.offset;
              if (offset == -1)
              {
                  offset_low  = -1;
          }
          if (file_version >= 8)
          {
 -            if (do_cpt_int(xd, "file_checksum_size", &(outputfiles[i].chksum_size),
 +            if (do_cpt_int(xd, "file_checksum_size", &outputfile.chksum_size,
                             list) != 0)
              {
                  return -1;
              }
 -            if (do_cpt_u_chars(xd, "file_checksum", 16, outputfiles[i].chksum, list) != 0)
 +            if (do_cpt_u_chars(xd, "file_checksum", 16, outputfile.chksum, list) != 0)
              {
                  return -1;
              }
          }
          else
          {
 -            outputfiles[i].chksum_size = -1;
 +            outputfile.chksum_size = -1;
          }
      }
      return 0;
  
  
  void write_checkpoint(const char *fn, gmx_bool bNumberAndKeep,
 -                      FILE *fplog, t_commrec *cr,
 +                      FILE *fplog, const t_commrec *cr,
                        ivec domdecCells, int nppnodes,
                        int eIntegrator, int simulation_part,
                        gmx_bool bExpanded, int elamstats,
                        t_state *state, ObservablesHistory *observablesHistory)
  {
      t_fileio            *fp;
 -    int                  file_version;
 -    char                *version;
 -    char                *btime;
 -    char                *buser;
 -    char                *bhost;
 -    int                  double_prec;
 -    char                *fprog;
      char                *fntemp; /* the temporary checkpoint file name */
 -    char                 timebuf[STRLEN];
      int                  npmenodes;
      char                 buf[1024], suffix[5+STEPSTRSIZE], sbuf[STEPSTRSIZE];
 -    gmx_file_position_t *outputfiles;
 -    int                  noutputfiles;
 -    char                *ftime;
      t_fileio            *ret;
  
      if (DOMAINDECOMP(cr))
      snew(fntemp, std::strlen(fn));
      std::strcpy(fntemp, fn);
  #endif
 +    char timebuf[STRLEN];
      gmx_format_current_time(timebuf, STRLEN);
  
      if (fplog)
      }
  
      /* Get offsets for open files */
 -    gmx_fio_get_output_file_positions(&outputfiles, &noutputfiles);
 +    auto outputfiles = gmx_fio_get_output_file_positions();
  
      fp = gmx_fio_open(fntemp, "w");
  
       * the same version, user, time, and build host!
       */
  
 -    version = gmx_strdup(gmx_version());
 -    btime   = gmx_strdup(BUILD_TIME);
 -    buser   = gmx_strdup(BUILD_USER);
 -    bhost   = gmx_strdup(BUILD_HOST);
 -
 -    double_prec = GMX_DOUBLE;
 -    fprog       = gmx_strdup(gmx::getProgramContext().fullBinaryPath());
 -
 -    ftime   = &(timebuf[0]);
 -
 -    int             nlambda     = (state->dfhist ? state->dfhist->nlambda : 0);
 -
 -    edsamhistory_t *edsamhist   = observablesHistory->edsamHistory.get();
 -    int             nED         = (edsamhist ? edsamhist->nED : 0);
 -
 -    swaphistory_t  *swaphist    = observablesHistory->swapHistory.get();
 -    int             eSwapCoords = (swaphist ? swaphist->eSwapCoords : eswapNO);
 -
 -    do_cpt_header(gmx_fio_getxdr(fp), FALSE, &file_version,
 -                  &version, &btime, &buser, &bhost, &double_prec, &fprog, &ftime,
 -                  &eIntegrator, &simulation_part, &step, &t, &nppnodes,
 -                  DOMAINDECOMP(cr) ? domdecCells : nullptr, &npmenodes,
 -                  &state->natoms, &state->ngtc, &state->nnhpres,
 -                  &state->nhchainlength, &nlambda, &state->flags, &flags_eks, &flags_enh, &flags_dfh, &flags_awhh,
 -                  &nED, &eSwapCoords,
 -                  nullptr);
 +    int                      nlambda     = (state->dfhist ? state->dfhist->nlambda : 0);
 +
 +    edsamhistory_t          *edsamhist   = observablesHistory->edsamHistory.get();
 +    int                      nED         = (edsamhist ? edsamhist->nED : 0);
 +
 +    swaphistory_t           *swaphist    = observablesHistory->swapHistory.get();
 +    int                      eSwapCoords = (swaphist ? swaphist->eSwapCoords : eswapNO);
 +
 +    CheckpointHeaderContents headerContents =
 +    {
 +        0, {0}, {0}, {0}, {0}, GMX_DOUBLE, {0}, {0},
 +        eIntegrator, simulation_part, step, t, nppnodes,
 +        {0}, npmenodes,
 +        state->natoms, state->ngtc, state->nnhpres,
 +        state->nhchainlength, nlambda, state->flags, flags_eks, flags_enh, flags_dfh, flags_awhh,
 +        nED, eSwapCoords
 +    };
 +    std::strcpy(headerContents.version, gmx_version());
 +    std::strcpy(headerContents.btime, BUILD_TIME);
 +    std::strcpy(headerContents.buser, BUILD_USER);
 +    std::strcpy(headerContents.bhost, BUILD_HOST);
 +    std::strcpy(headerContents.fprog, gmx::getProgramContext().fullBinaryPath());
 +    std::strcpy(headerContents.ftime, timebuf);
 +    if (DOMAINDECOMP(cr))
 +    {
 +        copy_ivec(domdecCells, headerContents.dd_nc);
 +    }
  
 -    sfree(version);
 -    sfree(btime);
 -    sfree(buser);
 -    sfree(bhost);
 -    sfree(fprog);
 +    do_cpt_header(gmx_fio_getxdr(fp), FALSE, nullptr, &headerContents);
  
      if ((do_cpt_state(gmx_fio_getxdr(fp), state->flags, state, nullptr) < 0)        ||
          (do_cpt_ekinstate(gmx_fio_getxdr(fp), flags_eks, &state->ekinstate, nullptr) < 0) ||
          (do_cpt_EDstate(gmx_fio_getxdr(fp), FALSE, nED, edsamhist, nullptr) < 0)      ||
          (do_cpt_awh(gmx_fio_getxdr(fp), FALSE, flags_awhh, state->awhHistory.get(), NULL) < 0) ||
          (do_cpt_swapstate(gmx_fio_getxdr(fp), FALSE, eSwapCoords, swaphist, nullptr) < 0) ||
 -        (do_cpt_files(gmx_fio_getxdr(fp), FALSE, &outputfiles, &noutputfiles, nullptr,
 -                      file_version) < 0))
 +        (do_cpt_files(gmx_fio_getxdr(fp), FALSE, &outputfiles, nullptr,
 +                      headerContents.file_version) < 0))
      {
          gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
      }
  
 -    do_cpt_footer(gmx_fio_getxdr(fp), file_version);
 +    do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
  
      /* we really, REALLY, want to make sure to physically write the checkpoint,
         and all the files it depends on, out to disk. Because we've
      }
  #endif  /* GMX_NO_RENAME */
  
 -    sfree(outputfiles);
      sfree(fntemp);
  
  #ifdef GMX_FAHCORE
  #endif /* end GMX_FAHCORE block */
  }
  
 -static void print_flag_mismatch(FILE *fplog, int sflags, int fflags)
 -{
 -    int i;
 -
 -    fprintf(fplog, "\nState entry mismatch between the simulation and the checkpoint file\n");
 -    fprintf(fplog, "Entries which are not present in the checkpoint file will not be updated\n");
 -    fprintf(fplog, "  %24s    %11s    %11s\n", "", "simulation", "checkpoint");
 -    for (i = 0; i < estNR; i++)
 -    {
 -        if ((sflags & (1<<i)) || (fflags & (1<<i)))
 -        {
 -            fprintf(fplog, "  %24s    %11s    %11s\n",
 -                    est_names[i],
 -                    (sflags & (1<<i)) ? "  present  " : "not present",
 -                    (fflags & (1<<i)) ? "  present  " : "not present");
 -        }
 -    }
 -}
 -
  static void check_int(FILE *fplog, const char *type, int p, int f, gmx_bool *mm)
  {
      FILE *fp = fplog ? fplog : stderr;
@@@ -2052,8 -2071,12 +2052,8 @@@ static void check_string(FILE *fplog, c
      }
  }
  
 -static void check_match(FILE *fplog,
 -                        char *version,
 -                        char *btime, char *buser, char *bhost, int double_prec,
 -                        char *fprog,
 -                        const t_commrec *cr, int npp_f, int npme_f,
 -                        const ivec dd_nc, const ivec dd_nc_f,
 +static void check_match(FILE *fplog, const t_commrec *cr, const ivec dd_nc,
 +                        const CheckpointHeaderContents &headerContents,
                          gmx_bool reproducibilityRequested)
  {
      /* Note that this check_string on the version will also print a message
       * message further down with reproducibilityRequested=TRUE.
       */
      gmx_bool versionDiffers = FALSE;
 -    check_string(fplog, "Version", gmx_version(), version, &versionDiffers);
 +    check_string(fplog, "Version", gmx_version(), headerContents.version, &versionDiffers);
  
      gmx_bool precisionDiffers = FALSE;
 -    check_int   (fplog, "Double prec.", GMX_DOUBLE, double_prec, &precisionDiffers);
 +    check_int   (fplog, "Double prec.", GMX_DOUBLE, headerContents.double_prec, &precisionDiffers);
      if (precisionDiffers)
      {
          const char msg_precision_difference[] =
  
      if (reproducibilityRequested)
      {
 -        check_string(fplog, "Build time", BUILD_TIME, btime, &mm);
 -        check_string(fplog, "Build user", BUILD_USER, buser, &mm);
 -        check_string(fplog, "Build host", BUILD_HOST, bhost, &mm);
 -        check_string(fplog, "Program name", gmx::getProgramContext().fullBinaryPath(), fprog, &mm);
 +        check_string(fplog, "Build time", BUILD_TIME, headerContents.btime, &mm);
 +        check_string(fplog, "Build user", BUILD_USER, headerContents.buser, &mm);
 +        check_string(fplog, "Build host", BUILD_HOST, headerContents.bhost, &mm);
 +        check_string(fplog, "Program name", gmx::getProgramContext().fullBinaryPath(), headerContents.fprog, &mm);
  
 -        check_int   (fplog, "#ranks", cr->nnodes, npp_f+npme_f, &mm);
 +        check_int   (fplog, "#ranks", cr->nnodes, headerContents.nnodes, &mm);
      }
  
      if (cr->nnodes > 1 && reproducibilityRequested)
      {
 -        check_int (fplog, "#PME-ranks", cr->npmenodes, npme_f, &mm);
 +        check_int (fplog, "#PME-ranks", cr->npmenodes, headerContents.npme, &mm);
  
          int npp = cr->nnodes;
          if (cr->npmenodes >= 0)
          {
              npp -= cr->npmenodes;
          }
 +        int npp_f = headerContents.nnodes;
 +        if (headerContents.npme >= 0)
 +        {
 +            npp_f -= headerContents.npme;
 +        }
          if (npp == npp_f)
          {
 -            check_int (fplog, "#DD-cells[x]", dd_nc[XX], dd_nc_f[XX], &mm);
 -            check_int (fplog, "#DD-cells[y]", dd_nc[YY], dd_nc_f[YY], &mm);
 -            check_int (fplog, "#DD-cells[z]", dd_nc[ZZ], dd_nc_f[ZZ], &mm);
 +            check_int (fplog, "#DD-cells[x]", dd_nc[XX], headerContents.dd_nc[XX], &mm);
 +            check_int (fplog, "#DD-cells[y]", dd_nc[YY], headerContents.dd_nc[YY], &mm);
 +            check_int (fplog, "#DD-cells[z]", dd_nc[ZZ], headerContents.dd_nc[ZZ], &mm);
          }
      }
  
          int        gmx_major;
          int        cpt_major;
          sscanf(gmx_version(), "%5d", &gmx_major);
 -        int        ret                 = sscanf(version, "%5d", &cpt_major);
 +        int        ret                 = sscanf(headerContents.version, "%5d", &cpt_major);
          gmx_bool   majorVersionDiffers = (ret < 1 || gmx_major != cpt_major);
  
          const char msg_major_version_difference[] =
  static void read_checkpoint(const char *fn, FILE **pfplog,
                              const t_commrec *cr,
                              const ivec dd_nc,
 -                            int eIntegrator, int *init_fep_state, gmx_int64_t *step, double *t,
 +                            int eIntegrator,
 +                            int *init_fep_state,
 +                            CheckpointHeaderContents *headerContents,
                              t_state *state, gmx_bool *bReadEkin,
                              ObservablesHistory *observablesHistory,
 -                            int *simulation_part,
                              gmx_bool bAppendOutputFiles, gmx_bool bForceAppend,
                              gmx_bool reproducibilityRequested)
  {
      t_fileio            *fp;
 -    int                  i, j, rc;
 -    int                  file_version;
 -    char                *version, *btime, *buser, *bhost, *fprog, *ftime;
 -    int                  double_prec;
 +    int                  j, rc;
      char                 buf[STEPSTRSIZE];
 -    int                  eIntegrator_f, nppnodes_f, npmenodes_f;
 -    ivec                 dd_nc_f;
 -    int                  natoms, ngtc, nnhpres, nhchainlength, nlambda, fflags, flags_eks, flags_enh, flags_dfh, flags_awhh;
 -    int                  nED, eSwapCoords;
      int                  ret;
 -    gmx_file_position_t *outputfiles;
 -    int                  nfiles;
      t_fileio            *chksum_file;
      FILE               * fplog = *pfplog;
      unsigned char        digest[16];
                                  dependent! */
  #endif
  
 -    const char *int_warn =
 -        "WARNING: The checkpoint file was generated with integrator %s,\n"
 -        "         while the simulation uses integrator %s\n\n";
 -
  #if !defined __native_client__ && !GMX_NATIVE_WINDOWS
      fl.l_type   = F_WRLCK;
      fl.l_whence = SEEK_SET;
  #endif
  
      fp = gmx_fio_open(fn, "r");
 -    do_cpt_header(gmx_fio_getxdr(fp), TRUE, &file_version,
 -                  &version, &btime, &buser, &bhost, &double_prec, &fprog, &ftime,
 -                  &eIntegrator_f, simulation_part, step, t,
 -                  &nppnodes_f, dd_nc_f, &npmenodes_f,
 -                  &natoms, &ngtc, &nnhpres, &nhchainlength, &nlambda,
 -                  &fflags, &flags_eks, &flags_enh, &flags_dfh, &flags_awhh,
 -                  &nED, &eSwapCoords, nullptr);
 +    do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, headerContents);
  
      if (bAppendOutputFiles &&
 -        file_version >= 13 && double_prec != GMX_DOUBLE)
 +        headerContents->file_version >= 13 && headerContents->double_prec != GMX_DOUBLE)
      {
          gmx_fatal(FARGS, "Output file appending requested, but the code and checkpoint file precision (single/double) don't match");
      }
      if (cr == nullptr || MASTER(cr))
      {
          fprintf(stderr, "\nReading checkpoint file %s generated: %s\n\n",
 -                fn, ftime);
 +                fn, headerContents->ftime);
      }
  
      /* This will not be written if we do appending, since fplog is still NULL then */
      {
          fprintf(fplog, "\n");
          fprintf(fplog, "Reading checkpoint file %s\n", fn);
 -        fprintf(fplog, "  file generated by:     %s\n", fprog);
 -        fprintf(fplog, "  file generated at:     %s\n", ftime);
 -        fprintf(fplog, "  GROMACS build time:    %s\n", btime);
 -        fprintf(fplog, "  GROMACS build user:    %s\n", buser);
 -        fprintf(fplog, "  GROMACS build host:    %s\n", bhost);
 -        fprintf(fplog, "  GROMACS double prec.:  %d\n", double_prec);
 -        fprintf(fplog, "  simulation part #:     %d\n", *simulation_part);
 -        fprintf(fplog, "  step:                  %s\n", gmx_step_str(*step, buf));
 -        fprintf(fplog, "  time:                  %f\n", *t);
 +        fprintf(fplog, "  file generated by:     %s\n", headerContents->fprog);
 +        fprintf(fplog, "  file generated at:     %s\n", headerContents->ftime);
 +        fprintf(fplog, "  GROMACS build time:    %s\n", headerContents->btime);
 +        fprintf(fplog, "  GROMACS build user:    %s\n", headerContents->buser);
 +        fprintf(fplog, "  GROMACS build host:    %s\n", headerContents->bhost);
 +        fprintf(fplog, "  GROMACS double prec.:  %d\n", headerContents->double_prec);
 +        fprintf(fplog, "  simulation part #:     %d\n", headerContents->simulation_part);
 +        fprintf(fplog, "  step:                  %s\n", gmx_step_str(headerContents->step, buf));
 +        fprintf(fplog, "  time:                  %f\n", headerContents->t);
          fprintf(fplog, "\n");
      }
  
 -    if (natoms != state->natoms)
 +    if (headerContents->natoms != state->natoms)
      {
 -        gmx_fatal(FARGS, "Checkpoint file is for a system of %d atoms, while the current system consists of %d atoms", natoms, state->natoms);
 +        gmx_fatal(FARGS, "Checkpoint file is for a system of %d atoms, while the current system consists of %d atoms", headerContents->natoms, state->natoms);
      }
 -    if (ngtc != state->ngtc)
 +    if (headerContents->ngtc != state->ngtc)
      {
 -        gmx_fatal(FARGS, "Checkpoint file is for a system of %d T-coupling groups, while the current system consists of %d T-coupling groups", ngtc, state->ngtc);
 +        gmx_fatal(FARGS, "Checkpoint file is for a system of %d T-coupling groups, while the current system consists of %d T-coupling groups", headerContents->ngtc, state->ngtc);
      }
 -    if (nnhpres != state->nnhpres)
 +    if (headerContents->nnhpres != state->nnhpres)
      {
 -        gmx_fatal(FARGS, "Checkpoint file is for a system of %d NH-pressure-coupling variables, while the current system consists of %d NH-pressure-coupling variables", nnhpres, state->nnhpres);
 +        gmx_fatal(FARGS, "Checkpoint file is for a system of %d NH-pressure-coupling variables, while the current system consists of %d NH-pressure-coupling variables", headerContents->nnhpres, state->nnhpres);
      }
  
      int nlambdaHistory = (state->dfhist ? state->dfhist->nlambda : 0);
 -    if (nlambda != nlambdaHistory)
 +    if (headerContents->nlambda != nlambdaHistory)
      {
 -        gmx_fatal(FARGS, "Checkpoint file is for a system with %d lambda states, while the current system consists of %d lambda states", nlambda, nlambdaHistory);
 +        gmx_fatal(FARGS, "Checkpoint file is for a system with %d lambda states, while the current system consists of %d lambda states", headerContents->nlambda, nlambdaHistory);
      }
  
 -    init_gtc_state(state, state->ngtc, state->nnhpres, nhchainlength); /* need to keep this here to keep the tpr format working */
 +    init_gtc_state(state, state->ngtc, state->nnhpres, headerContents->nhchainlength); /* need to keep this here to keep the tpr format working */
      /* write over whatever was read; we use the number of Nose-Hoover chains from the checkpoint */
  
 -    if (eIntegrator_f != eIntegrator)
 +    if (headerContents->eIntegrator != eIntegrator)
      {
 -        if (MASTER(cr))
 -        {
 -            fprintf(stderr, int_warn, EI(eIntegrator_f), EI(eIntegrator));
 -        }
 -        if (bAppendOutputFiles)
 -        {
 -            gmx_fatal(FARGS,
 -                      "Output file appending requested, but input/checkpoint integrators do not match.\n"
 -                      "Stopping the run to prevent you from ruining all your data...\n"
 -                      "If you _really_ know what you are doing, try with the -noappend option.\n");
 -        }
 -        if (fplog)
 -        {
 -            fprintf(fplog, int_warn, EI(eIntegrator_f), EI(eIntegrator));
 -        }
 +        gmx_fatal(FARGS, "Cannot change integrator during a checkpoint restart. Perhaps you should make a new .tpr with grompp -f new.mdp -t %s", fn);
      }
  
 -    if (fflags != state->flags)
 +    if (headerContents->flags_state != state->flags)
      {
 -
 -        if (MASTER(cr))
 -        {
 -            if (bAppendOutputFiles)
 -            {
 -                gmx_fatal(FARGS,
 -                          "Output file appending requested, but input and checkpoint states are not identical.\n"
 -                          "Stopping the run to prevent you from ruining all your data...\n"
 -                          "You can try with the -noappend option, and get more info in the log file.\n");
 -            }
 -
 -            if (getenv("GMX_ALLOW_CPT_MISMATCH") == nullptr)
 -            {
 -                gmx_fatal(FARGS, "You seem to have switched ensemble, integrator, T and/or P-coupling algorithm between the cpt and tpr file, or switched GROMACS version. The recommended way of doing this is passing the cpt file to grompp (with option -t) instead of to mdrun. If you know what you are doing, you can override this error by setting the env.var. GMX_ALLOW_CPT_MISMATCH");
 -            }
 -            else
 -            {
 -                fprintf(stderr,
 -                        "WARNING: The checkpoint state entries do not match the simulation,\n"
 -                        "         see the log file for details\n\n");
 -            }
 -        }
 -
 -        if (fplog)
 -        {
 -            print_flag_mismatch(fplog, state->flags, fflags);
 -        }
 +        gmx_fatal(FARGS, "Cannot change a simulation algorithm during a checkpoint restart. Perhaps you should make a new .tpr with grompp -f new.mdp -t %s", fn);
      }
 -    else
 +
 +    if (MASTER(cr))
      {
 -        if (MASTER(cr))
 -        {
 -            check_match(fplog, version, btime, buser, bhost, double_prec, fprog,
 -                        cr, nppnodes_f, npmenodes_f, dd_nc, dd_nc_f,
 -                        reproducibilityRequested);
 -        }
 +        check_match(fplog, cr, dd_nc, *headerContents,
 +                    reproducibilityRequested);
      }
 -    ret             = do_cpt_state(gmx_fio_getxdr(fp), fflags, state, nullptr);
 +
 +    ret             = do_cpt_state(gmx_fio_getxdr(fp), headerContents->flags_state, state, nullptr);
      *init_fep_state = state->fep_state;  /* there should be a better way to do this than setting it here.
                                              Investigate for 5.0. */
      if (ret)
      {
          cp_error();
      }
 -    ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), flags_eks, &state->ekinstate, nullptr);
 +    ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents->flags_eks, &state->ekinstate, nullptr);
      if (ret)
      {
          cp_error();
      }
 -    *bReadEkin = ((flags_eks & (1<<eeksEKINH)) || (flags_eks & (1<<eeksEKINF)) || (flags_eks & (1<<eeksEKINO)) ||
 -                  ((flags_eks & (1<<eeksEKINSCALEF)) | (flags_eks & (1<<eeksEKINSCALEH)) | (flags_eks & (1<<eeksVSCALE))));
 +    *bReadEkin = ((headerContents->flags_eks & (1<<eeksEKINH)) ||
 +                  (headerContents->flags_eks & (1<<eeksEKINF)) ||
 +                  (headerContents->flags_eks & (1<<eeksEKINO)) ||
 +                  ((headerContents->flags_eks & (1<<eeksEKINSCALEF)) |
 +                   (headerContents->flags_eks & (1<<eeksEKINSCALEH)) |
 +                   (headerContents->flags_eks & (1<<eeksVSCALE))));
  
 -    if (flags_enh && observablesHistory->energyHistory == nullptr)
 +    if (headerContents->flags_enh && observablesHistory->energyHistory == nullptr)
      {
          observablesHistory->energyHistory = std::unique_ptr<energyhistory_t>(new energyhistory_t {});
      }
      ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE,
 -                          flags_enh, observablesHistory->energyHistory.get(), nullptr);
 +                          headerContents->flags_enh, observablesHistory->energyHistory.get(), nullptr);
      if (ret)
      {
          cp_error();
      }
  
 -    if (file_version < 6)
 +    if (headerContents->file_version < 6)
      {
          const char *warn = "Reading checkpoint file in old format, assuming that the run that generated this file started at step 0, if this is not the case the averages stored in the energy file will be incorrect.";
  
          {
              fprintf(fplog, "\nWARNING: %s\n\n", warn);
          }
 -        observablesHistory->energyHistory->nsum     = *step;
 -        observablesHistory->energyHistory->nsum_sim = *step;
 +        observablesHistory->energyHistory->nsum     = headerContents->step;
 +        observablesHistory->energyHistory->nsum_sim = headerContents->step;
      }
  
 -    ret = do_cpt_df_hist(gmx_fio_getxdr(fp), flags_dfh, nlambda, &state->dfhist, nullptr);
 +    ret = do_cpt_df_hist(gmx_fio_getxdr(fp), headerContents->flags_dfh, headerContents->nlambda, &state->dfhist, nullptr);
      if (ret)
      {
          cp_error();
      }
  
 -    if (nED > 0 && observablesHistory->edsamHistory == nullptr)
 +    if (headerContents->nED > 0 && observablesHistory->edsamHistory == nullptr)
      {
          observablesHistory->edsamHistory = std::unique_ptr<edsamhistory_t>(new edsamhistory_t {});
      }
 -    ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, nED, observablesHistory->edsamHistory.get(), nullptr);
 +    ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, headerContents->nED, observablesHistory->edsamHistory.get(), nullptr);
      if (ret)
      {
          cp_error();
      }
  
 -    if (flags_awhh != 0 && state->awhHistory == nullptr)
 +    if (headerContents->flags_awhh != 0 && state->awhHistory == nullptr)
      {
          state->awhHistory = std::shared_ptr<gmx::AwhHistory>(new gmx::AwhHistory());
      }
      ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE,
 -                     flags_awhh, state->awhHistory.get(), NULL);
 +                     headerContents->flags_awhh, state->awhHistory.get(), NULL);
      if (ret)
      {
          cp_error();
      }
  
 -    if (eSwapCoords != eswapNO && observablesHistory->swapHistory == nullptr)
 +    if (headerContents->eSwapCoords != eswapNO && observablesHistory->swapHistory == nullptr)
      {
          observablesHistory->swapHistory = std::unique_ptr<swaphistory_t>(new swaphistory_t {});
      }
 -    ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, eSwapCoords, observablesHistory->swapHistory.get(), nullptr);
 +    ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents->eSwapCoords, observablesHistory->swapHistory.get(), nullptr);
      if (ret)
      {
          cp_error();
      }
  
 -    ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, &nfiles, nullptr, file_version);
 +    std::vector<gmx_file_position_t> outputfiles;
 +    ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, nullptr, headerContents->file_version);
      if (ret)
      {
          cp_error();
      }
  
 -    ret = do_cpt_footer(gmx_fio_getxdr(fp), file_version);
 +    ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents->file_version);
      if (ret)
      {
          cp_error();
          gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
      }
  
 -    sfree(fprog);
 -    sfree(ftime);
 -    sfree(btime);
 -    sfree(buser);
 -    sfree(bhost);
 -
      /* If the user wants to append to output files,
       * we use the file pointer positions of the output files stored
       * in the checkpoint file and truncate the files such that any frames
       */
      if (bAppendOutputFiles)
      {
 +        if (outputfiles.empty())
 +        {
 +            gmx_fatal(FARGS, "No names of output files were recorded in the checkpoint");
 +        }
          if (fn2ftp(outputfiles[0].filename) != efLOG)
          {
              /* make sure first file is log file so that it is OK to use it for
              gmx_fatal(FARGS, "The first output file should always be the log "
                        "file but instead is: %s. Cannot do appending because of this condition.", outputfiles[0].filename);
          }
 -        for (i = 0; i < nfiles; i++)
 +        bool firstFile = true;
 +        for (const auto &outputfile : outputfiles)
          {
 -            if (outputfiles[i].offset < 0)
 +            if (outputfile.offset < 0)
              {
                  gmx_fatal(FARGS, "The original run wrote a file called '%s' which "
                            "is larger than 2 GB, but mdrun did not support large file"
                            " offsets. Can not append. Run mdrun with -noappend",
 -                          outputfiles[i].filename);
 +                          outputfile.filename);
              }
  #ifdef GMX_FAHCORE
 -            chksum_file = gmx_fio_open(outputfiles[i].filename, "a");
 +            chksum_file = gmx_fio_open(outputfile.filename, "a");
  
  #else
 -            chksum_file = gmx_fio_open(outputfiles[i].filename, "r+");
 +            chksum_file = gmx_fio_open(outputfile.filename, "r+");
  
              /* lock log file */
 -            if (i == 0)
 +            if (firstFile)
              {
                  /* Note that there are systems where the lock operation
                   * will succeed, but a second process can also lock the file.
                          }
                          else
                          {
 -                            fprintf(stderr, "\nNOTE: File locking is not supported on this system, will not lock %s\n\n", outputfiles[i].filename);
 +                            fprintf(stderr, "\nNOTE: File locking is not supported on this system, will not lock %s\n\n", outputfile.filename);
                              if (fplog)
                              {
 -                                fprintf(fplog, "\nNOTE: File locking not supported on this system, will not lock %s\n\n", outputfiles[i].filename);
 +                                fprintf(fplog, "\nNOTE: File locking not supported on this system, will not lock %s\n\n", outputfile.filename);
                              }
                          }
                      }
                      else if (errno == EACCES || errno == EAGAIN)
                      {
                          gmx_fatal(FARGS, "Failed to lock: %s. Already running "
 -                                  "simulation?", outputfiles[i].filename);
 +                                  "simulation?", outputfile.filename);
                      }
                      else
                      {
                          gmx_fatal(FARGS, "Failed to lock: %s. %s.",
 -                                  outputfiles[i].filename, std::strerror(errno));
 +                                  outputfile.filename, std::strerror(errno));
                      }
                  }
              }
  
              /* compute md5 chksum */
 -            if (outputfiles[i].chksum_size != -1)
 +            if (outputfile.chksum_size != -1)
              {
 -                if (gmx_fio_get_file_md5(chksum_file, outputfiles[i].offset,
 -                                         digest) != outputfiles[i].chksum_size) /*at the end of the call the file position is at the end of the file*/
 +                if (gmx_fio_get_file_md5(chksum_file, outputfile.offset,
 +                                         digest) != outputfile.chksum_size) /*at the end of the call the file position is at the end of the file*/
                  {
                      gmx_fatal(FARGS, "Can't read %d bytes of '%s' to compute checksum. The file has been replaced or its contents have been modified. Cannot do appending because of this condition.",
 -                              outputfiles[i].chksum_size,
 -                              outputfiles[i].filename);
 +                              outputfile.chksum_size,
 +                              outputfile.filename);
                  }
              }
 -            if (i == 0)  /*log file needs to be seeked in case we need to truncate (other files are truncated below)*/
 +            /* log file needs to be seeked in case we need to truncate
 +             * (other files are truncated below)*/
 +            if (firstFile)
              {
 -                if (gmx_fio_seek(chksum_file, outputfiles[i].offset))
 +                if (gmx_fio_seek(chksum_file, outputfile.offset))
                  {
                      gmx_fatal(FARGS, "Seek error! Failed to truncate log-file: %s.", std::strerror(errno));
                  }
              }
  #endif
  
 -            if (i == 0) /*open log file here - so that lock is never lifted
 -                           after chksum is calculated */
 +            /* open log file here - so that lock is never lifted
 +             * after chksum is calculated */
 +            if (firstFile)
              {
                  *pfplog = gmx_fio_getfp(chksum_file);
              }
              }
  #ifndef GMX_FAHCORE
              /* compare md5 chksum */
 -            if (outputfiles[i].chksum_size != -1 &&
 -                memcmp(digest, outputfiles[i].chksum, 16) != 0)
 +            if (outputfile.chksum_size != -1 &&
 +                memcmp(digest, outputfile.chksum, 16) != 0)
              {
                  if (debug)
                  {
 -                    fprintf(debug, "chksum for %s: ", outputfiles[i].filename);
 +                    fprintf(debug, "chksum for %s: ", outputfile.filename);
                      for (j = 0; j < 16; j++)
                      {
                          fprintf(debug, "%02x", digest[j]);
                      fprintf(debug, "\n");
                  }
                  gmx_fatal(FARGS, "Checksum wrong for '%s'. The file has been replaced or its contents have been modified. Cannot do appending because of this condition.",
 -                          outputfiles[i].filename);
 +                          outputfile.filename);
              }
  #endif
  
  
 -            if (i != 0) /*log file is already seeked to correct position */
 +            if (!firstFile) /* log file is already seeked to correct position */
              {
  #if !GMX_NATIVE_WINDOWS || !defined(GMX_FAHCORE)
                  /* For FAHCORE, we do this elsewhere*/
 -                rc = gmx_truncate(outputfiles[i].filename, outputfiles[i].offset);
 +                rc = gmx_truncate(outputfile.filename, outputfile.offset);
                  if (rc != 0)
                  {
 -                    gmx_fatal(FARGS, "Truncation of file %s failed. Cannot do appending because of this failure.", outputfiles[i].filename);
 +                    gmx_fatal(FARGS, "Truncation of file %s failed. Cannot do appending because of this failure.", outputfile.filename);
                  }
  #endif
              }
 +            firstFile = false;
          }
      }
 -
 -    sfree(outputfiles);
  }
  
  
@@@ -2517,29 -2589,30 +2517,29 @@@ void load_checkpoint(const char *fn, FI
                       gmx_bool bAppend, gmx_bool bForceAppend,
                       gmx_bool reproducibilityRequested)
  {
 -    gmx_int64_t     step;
 -    double          t;
 -
 +    CheckpointHeaderContents headerContents;
      if (SIMMASTER(cr))
      {
          /* Read the state from the checkpoint file */
          read_checkpoint(fn, fplog,
                          cr, dd_nc,
 -                        ir->eI, &(ir->fepvals->init_fep_state), &step, &t,
 +                        ir->eI, &(ir->fepvals->init_fep_state),
 +                        &headerContents,
                          state, bReadEkin, observablesHistory,
 -                        &ir->simulation_part, bAppend, bForceAppend,
 +                        bAppend, bForceAppend,
                          reproducibilityRequested);
      }
      if (PAR(cr))
      {
 -        gmx_bcast(sizeof(step), &step, cr);
 +        gmx_bcast(sizeof(headerContents.step), &headerContents.step, cr);
          gmx_bcast(sizeof(*bReadEkin), bReadEkin, cr);
      }
      ir->bContinuation    = TRUE;
      if (ir->nsteps >= 0)
      {
 -        ir->nsteps          += ir->init_step - step;
 +        ir->nsteps          += ir->init_step - headerContents.step;
      }
 -    ir->init_step        = step;
 +    ir->init_step        = headerContents.step;
      ir->simulation_part += 1;
  }
  
@@@ -2547,6 -2620,17 +2547,6 @@@ void read_checkpoint_part_and_step(cons
                                     int         *simulation_part,
                                     gmx_int64_t *step)
  {
 -    int       file_version;
 -    char     *version, *btime, *buser, *bhost, *fprog, *ftime;
 -    int       double_prec;
 -    int       eIntegrator;
 -    int       nppnodes, npme;
 -    ivec      dd_nc;
 -    int       nlambda;
 -    int       flags_eks, flags_enh, flags_dfh, flags_awhh;
 -    double    t;
 -    t_state   state;
 -    int       nED, eSwapCoords;
      t_fileio *fp;
  
      if (filename == nullptr ||
          return;
      }
  
 -    /* Not calling initializing state before use is nasty, but all we
 -       do is read into its member variables and throw the struct away
 -       again immediately. */
 -
 -    do_cpt_header(gmx_fio_getxdr(fp), TRUE, &file_version,
 -                  &version, &btime, &buser, &bhost, &double_prec, &fprog, &ftime,
 -                  &eIntegrator, simulation_part, step, &t, &nppnodes, dd_nc, &npme,
 -                  &state.natoms, &state.ngtc, &state.nnhpres, &state.nhchainlength,
 -                  &nlambda, &state.flags, &flags_eks, &flags_enh, &flags_dfh, &flags_awhh,
 -                  &nED, &eSwapCoords, nullptr);
 -
 +    CheckpointHeaderContents headerContents;
 +    do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, &headerContents);
      gmx_fio_close(fp);
 +    *simulation_part = headerContents.simulation_part;
 +    *step            = headerContents.step;
  }
  
  static void read_checkpoint_data(t_fileio *fp, int *simulation_part,
                                   gmx_int64_t *step, double *t, t_state *state,
 -                                 int *nfiles, gmx_file_position_t **outputfiles)
 -{
 -    int                  file_version;
 -    char                *version, *btime, *buser, *bhost, *fprog, *ftime;
 -    int                  double_prec;
 -    int                  eIntegrator;
 -    int                  nppnodes, npme;
 -    ivec                 dd_nc;
 -    int                  nlambda;
 -    int                  flags_eks, flags_enh, flags_dfh, flags_awhh;
 -    int                  nED, eSwapCoords;
 -    int                  nfiles_loc;
 -    gmx_file_position_t *files_loc = nullptr;
 -    int                  ret;
 -
 -    do_cpt_header(gmx_fio_getxdr(fp), TRUE, &file_version,
 -                  &version, &btime, &buser, &bhost, &double_prec, &fprog, &ftime,
 -                  &eIntegrator, simulation_part, step, t, &nppnodes, dd_nc, &npme,
 -                  &state->natoms, &state->ngtc, &state->nnhpres, &state->nhchainlength,
 -                  &nlambda, &state->flags, &flags_eks, &flags_enh, &flags_dfh, &flags_awhh,
 -                  &nED, &eSwapCoords, nullptr);
 -    ret =
 +                                 std::vector<gmx_file_position_t> *outputfiles)
 +{
 +    int                      ret;
 +
 +    CheckpointHeaderContents headerContents;
 +    do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, &headerContents);
 +    *simulation_part     = headerContents.simulation_part;
 +    *step                = headerContents.step;
 +    *t                   = headerContents.t;
 +    state->natoms        = headerContents.natoms;
 +    state->ngtc          = headerContents.ngtc;
 +    state->nnhpres       = headerContents.nnhpres;
 +    state->nhchainlength = headerContents.nhchainlength;
 +    state->flags         = headerContents.flags_state;
 +    ret                  =
          do_cpt_state(gmx_fio_getxdr(fp), state->flags, state, nullptr);
      if (ret)
      {
          cp_error();
      }
 -    ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), flags_eks, &state->ekinstate, nullptr);
 +    ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents.flags_eks, &state->ekinstate, nullptr);
      if (ret)
      {
          cp_error();
  
      energyhistory_t enerhist;
      ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE,
 -                          flags_enh, &enerhist, nullptr);
 +                          headerContents.flags_enh, &enerhist, nullptr);
      if (ret)
      {
          cp_error();
      }
 -    ret = do_cpt_df_hist(gmx_fio_getxdr(fp), flags_dfh, nlambda, &state->dfhist, nullptr);
 +    ret = do_cpt_df_hist(gmx_fio_getxdr(fp), headerContents.flags_dfh, headerContents.nlambda, &state->dfhist, nullptr);
      if (ret)
      {
          cp_error();
      }
  
      edsamhistory_t edsamhist = {};
 -    ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, nED, &edsamhist, nullptr);
 +    ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, headerContents.nED, &edsamhist, nullptr);
      if (ret)
      {
          cp_error();
      }
  
      ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE,
 -                     flags_awhh, state->awhHistory.get(), NULL);
 +                     headerContents.flags_awhh, state->awhHistory.get(), NULL);
      if (ret)
      {
          cp_error();
      }
  
      swaphistory_t swaphist = {};
 -    ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, eSwapCoords, &swaphist, nullptr);
 +    ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents.eSwapCoords, &swaphist, nullptr);
      if (ret)
      {
          cp_error();
      }
  
      ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE,
 -                       &files_loc,
 -                       &nfiles_loc,
 -                       nullptr, file_version);
 -    if (outputfiles != nullptr)
 -    {
 -        *outputfiles = files_loc;
 -    }
 -    else
 -    {
 -        sfree(files_loc);
 -    }
 -    if (nfiles != nullptr)
 -    {
 -        *nfiles = nfiles_loc;
 -    }
 +                       outputfiles,
 +                       nullptr, headerContents.file_version);
  
      if (ret)
      {
          cp_error();
      }
  
 -    ret = do_cpt_footer(gmx_fio_getxdr(fp), file_version);
 +    ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
      if (ret)
      {
          cp_error();
      }
 -
 -    sfree(fprog);
 -    sfree(ftime);
 -    sfree(btime);
 -    sfree(buser);
 -    sfree(bhost);
 -}
 -
 -void
 -read_checkpoint_state(const char *fn, int *simulation_part,
 -                      gmx_int64_t *step, double *t, t_state *state)
 -{
 -    t_fileio *fp;
 -
 -    fp = gmx_fio_open(fn, "r");
 -    read_checkpoint_data(fp, simulation_part, step, t, state, nullptr, nullptr);
 -    if (gmx_fio_close(fp) != 0)
 -    {
 -        gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
 -    }
  }
  
  void read_checkpoint_trxframe(t_fileio *fp, t_trxframe *fr)
  {
 -    t_state         state;
 -    int             simulation_part;
 -    gmx_int64_t     step;
 -    double          t;
 +    t_state                          state;
 +    int                              simulation_part;
 +    gmx_int64_t                      step;
 +    double                           t;
  
 -    read_checkpoint_data(fp, &simulation_part, &step, &t, &state, nullptr, nullptr);
 +    std::vector<gmx_file_position_t> outputfiles;
 +    read_checkpoint_data(fp, &simulation_part, &step, &t, &state, &outputfiles);
  
      fr->natoms  = state.natoms;
      fr->bStep   = TRUE;
  void list_checkpoint(const char *fn, FILE *out)
  {
      t_fileio            *fp;
 -    int                  file_version;
 -    char                *version, *btime, *buser, *bhost, *fprog, *ftime;
 -    int                  double_prec;
 -    int                  eIntegrator, simulation_part, nppnodes, npme;
 -    gmx_int64_t          step;
 -    double               t;
 -    ivec                 dd_nc;
 -    int                  nlambda;
 -    int                  flags_eks, flags_enh, flags_dfh, flags_awhh;;
 -    int                  nED, eSwapCoords;
      int                  ret;
 -    gmx_file_position_t *outputfiles;
 -    int                  nfiles;
  
      t_state              state;
  
      fp = gmx_fio_open(fn, "r");
 -    do_cpt_header(gmx_fio_getxdr(fp), TRUE, &file_version,
 -                  &version, &btime, &buser, &bhost, &double_prec, &fprog, &ftime,
 -                  &eIntegrator, &simulation_part, &step, &t, &nppnodes, dd_nc, &npme,
 -                  &state.natoms, &state.ngtc, &state.nnhpres, &state.nhchainlength,
 -                  &nlambda, &state.flags,
 -                  &flags_eks, &flags_enh, &flags_dfh, &flags_awhh, &nED, &eSwapCoords,
 -                  out);
 -    ret = do_cpt_state(gmx_fio_getxdr(fp), state.flags, &state, out);
 +    CheckpointHeaderContents headerContents;
 +    do_cpt_header(gmx_fio_getxdr(fp), TRUE, out, &headerContents);
 +    state.natoms        = headerContents.natoms;
 +    state.ngtc          = headerContents.ngtc;
 +    state.nnhpres       = headerContents.nnhpres;
 +    state.nhchainlength = headerContents.nhchainlength;
 +    state.flags         = headerContents.flags_state;
 +    ret                 = do_cpt_state(gmx_fio_getxdr(fp), state.flags, &state, out);
      if (ret)
      {
          cp_error();
      }
 -    ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), flags_eks, &state.ekinstate, out);
 +    ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents.flags_eks, &state.ekinstate, out);
      if (ret)
      {
          cp_error();
  
      energyhistory_t enerhist;
      ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE,
 -                          flags_enh, &enerhist, out);
 +                          headerContents.flags_enh, &enerhist, out);
  
      if (ret == 0)
      {
          ret = do_cpt_df_hist(gmx_fio_getxdr(fp),
 -                             flags_dfh, nlambda, &state.dfhist, out);
 +                             headerContents.flags_dfh, headerContents.nlambda, &state.dfhist, out);
      }
  
      if (ret == 0)
      {
          edsamhistory_t edsamhist = {};
 -        ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, nED, &edsamhist, out);
 +        ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, headerContents.nED, &edsamhist, out);
      }
  
      if (ret == 0)
      {
          ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE,
 -                         flags_awhh, state.awhHistory.get(), out);
 +                         headerContents.flags_awhh, state.awhHistory.get(), out);
      }
  
      if (ret == 0)
      {
          swaphistory_t swaphist = {};
 -        ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, eSwapCoords, &swaphist, out);
 +        ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents.eSwapCoords, &swaphist, out);
      }
  
      if (ret == 0)
      {
 -        ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, &nfiles, out, file_version);
 +        std::vector<gmx_file_position_t> outputfiles;
 +        ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, out, headerContents.file_version);
      }
  
      if (ret == 0)
      {
 -        ret = do_cpt_footer(gmx_fio_getxdr(fp), file_version);
 +        ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
      }
  
      if (ret)
  
  /* This routine cannot print tons of data, since it is called before the log file is opened. */
  void
 -read_checkpoint_simulation_part_and_filenames(t_fileio             *fp,
 -                                              int                  *simulation_part,
 -                                              int                  *nfiles,
 -                                              gmx_file_position_t **outputfiles)
 +read_checkpoint_simulation_part_and_filenames(t_fileio                         *fp,
 +                                              int                              *simulation_part,
 +                                              std::vector<gmx_file_position_t> *outputfiles)
  {
      gmx_int64_t step = 0;
      double      t;
      t_state     state;
  
 -    read_checkpoint_data(fp, simulation_part, &step, &t, &state,
 -                         nfiles, outputfiles);
 +    read_checkpoint_data(fp, simulation_part, &step, &t, &state, outputfiles);
      if (gmx_fio_close(fp) != 0)
      {
          gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
index 44c5c6e1fc8286598654ca6b421bc8f3821b9a75,2002f94ca006b3734b4c4eda5f836cfef9f28b04..503dad983e35b3cae1177ee4157e876604c0a697
@@@ -55,7 -55,6 +55,7 @@@
  #include "gromacs/utility/exceptions.h"
  #include "gromacs/utility/fatalerror.h"
  #include "gromacs/utility/futil.h"
 +#include "gromacs/utility/gmxassert.h"
  #include "gromacs/utility/gmxmpi.h"
  #include "gromacs/utility/path.h"
  #include "gromacs/utility/programcontext.h"
@@@ -239,53 -238,56 +239,61 @@@ void gmx_log_close(FILE *fp
  }
  
  // TODO move this to multi-sim module
 -void init_multisystem(t_commrec *cr, int nsim, char **multidirs,
 -                      int nfile, const t_filenm fnm[])
 +gmx_multisim_t *init_multisystem(MPI_Comm                         comm,
 +                                 gmx::ArrayRef<const std::string> multidirs)
  {
      gmx_multisim_t *ms;
 -    int             nnodes, nnodpersim, sim, i;
  #if GMX_MPI
      MPI_Group       mpi_group_world;
      int            *rank;
  #endif
  
-     if (multidirs.size() <= 1)
 -#if !GMX_MPI
 -    if (nsim > 1)
++    if (multidirs.empty())
      {
 -        gmx_fatal(FARGS, "This binary is compiled without MPI support, can not do multiple simulations.");
 +        return nullptr;
 +    }
-     if (!GMX_LIB_MPI && multidirs.size() > 1)
++
++    if (!GMX_LIB_MPI && multidirs.size() >= 1)
 +    {
 +        gmx_fatal(FARGS, "mdrun -multidir is only supported when GROMACS has been "
 +                  "configured with a proper external MPI library.");
      }
 -#endif
  
 -    if (nsim == 1)
++    if (multidirs.size() == 1)
+     {
+         /* NOTE: It would be nice if this special case worked, but this requires checks/tests. */
 -        gmx_fatal(FARGS, "To run mdrun in multiple simulation mode, more then one actual simulation is required. The single simulation case is not supported.");
++        gmx_fatal(FARGS, "To run mdrun in multiple simulation mode, more then one "
++                  "actual simulation is required. The single simulation case is not supported.");
+     }
 -    nnodes  = cr->nnodes;
 -    if (nnodes % nsim != 0)
 +#if GMX_MPI
 +    int numRanks;
 +    MPI_Comm_size(comm, &numRanks);
 +    if (numRanks % multidirs.size() != 0)
      {
 -        gmx_fatal(FARGS, "The number of ranks (%d) is not a multiple of the number of simulations (%d)", nnodes, nsim);
 +        gmx_fatal(FARGS, "The number of ranks (%d) is not a multiple of the number of simulations (%zu)", numRanks, multidirs.size());
      }
  
 -    nnodpersim = nnodes/nsim;
 -    sim        = cr->nodeid/nnodpersim;
 +    int numRanksPerSim = numRanks/multidirs.size();
 +    int rankWithinComm;
 +    MPI_Comm_rank(comm, &rankWithinComm);
  
      if (debug)
      {
 -        fprintf(debug, "We have %d simulations, %d ranks per simulation, local simulation is %d\n", nsim, nnodpersim, sim);
 +        fprintf(debug, "We have %zu simulations, %d ranks per simulation, local simulation is %d\n", multidirs.size(), numRanksPerSim, rankWithinComm/numRanksPerSim);
      }
  
 -    snew(ms, 1);
 -    cr->ms   = ms;
 -    ms->nsim = nsim;
 -    ms->sim  = sim;
 -#if GMX_MPI
 +    ms       = new gmx_multisim_t;
 +    ms->nsim = multidirs.size();
 +    ms->sim  = rankWithinComm/numRanksPerSim;
      /* Create a communicator for the master nodes */
      snew(rank, ms->nsim);
 -    for (i = 0; i < ms->nsim; i++)
 +    for (int i = 0; i < ms->nsim; i++)
      {
 -        rank[i] = i*nnodpersim;
 +        rank[i] = i*numRanksPerSim;
      }
 -    MPI_Comm_group(MPI_COMM_WORLD, &mpi_group_world);
 -    MPI_Group_incl(mpi_group_world, nsim, rank, &ms->mpi_group_masters);
 +    MPI_Comm_group(comm, &mpi_group_world);
 +    MPI_Group_incl(mpi_group_world, ms->nsim, rank, &ms->mpi_group_masters);
      sfree(rank);
      MPI_Comm_create(MPI_COMM_WORLD, ms->mpi_group_masters,
                      &ms->mpi_comm_masters);
      ms->mpb->dbuf_alloc  = 0;
  #endif
  
 +    // TODO This should throw upon error
 +    gmx_chdir(multidirs[ms->sim].c_str());
 +#else
 +    GMX_UNUSED_VALUE(comm);
 +    ms = nullptr;
  #endif
  
 -    /* Reduce the intra-simulation communication */
 -    cr->sim_nodeid = cr->nodeid % nnodpersim;
 -    cr->nnodes     = nnodpersim;
 -#if GMX_MPI
 -    MPI_Comm_split(MPI_COMM_WORLD, sim, cr->sim_nodeid, &cr->mpi_comm_mysim);
 -    cr->mpi_comm_mygroup = cr->mpi_comm_mysim;
 -    cr->nodeid           = cr->sim_nodeid;
 -#endif
 -
 -    if (debug)
 -    {
 -        fprintf(debug, "This is simulation %d", cr->ms->sim);
 -        if (PAR(cr))
 -        {
 -            fprintf(debug, ", local number of ranks %d, local rank ID %d",
 -                    cr->nnodes, cr->sim_nodeid);
 -        }
 -        fprintf(debug, "\n\n");
 -    }
 +    return ms;
 +}
  
 -    if (multidirs)
 -    {
 -        if (debug)
 -        {
 -            fprintf(debug, "Changing to directory %s\n", multidirs[cr->ms->sim]);
 -        }
 -        gmx_chdir(multidirs[cr->ms->sim]);
 -    }
 -    else
 +void done_multisim(gmx_multisim_t *ms)
 +{
 +    if (nullptr != ms)
      {
 -        try
 -        {
 -            std::string rankString = gmx::formatString("%d", cr->ms->sim);
 -            /* Patch output and tpx, cpt and rerun input file names */
 -            for (i = 0; (i < nfile); i++)
 -            {
 -                /* Because of possible multiple extensions per type we must look
 -                 * at the actual file name for rerun. */
 -                if (is_output(&fnm[i]) ||
 -                    fnm[i].ftp == efTPR || fnm[i].ftp == efCPT ||
 -                    strcmp(fnm[i].opt, "-rerun") == 0)
 -                {
 -                    std::string newFileName = gmx::Path::concatenateBeforeExtension(fnm[i].fns[0], rankString);
 -                    sfree(fnm[i].fns[0]);
 -                    fnm[i].fns[0] = gmx_strdup(newFileName.c_str());
 -                }
 -            }
 -        }
 -        catch (gmx::GromacsException &e)
 -        {
 -            e.prependContext(gmx::formatString("Failed to modify mdrun -multi filename to add per-simulation suffix. You could perhaps reorganize your files and try mdrun -multidir.\n"));
 -            throw;
 -        }
 +        done_mpi_in_place_buf(ms->mpb);
 +        delete ms;
      }
  }
index c24d9018db46eb9cfde4f065d4b90756cb699e82,caff4df72c31c6f2e6497ad64220266c560d669a..9806b4f9c531f36efdb2b5747a005c2b696dd3d2
@@@ -938,10 -938,13 +938,13 @@@ static double get_pull_coord_deviation(
  
      pcrd = &pull->coord[coord_ind];
  
-     get_pull_coord_distance(pull, coord_ind, pbc);
+     /* Update the reference value before computing the distance,
+      * since it is used in the distance computation with periodic pulling.
+      */
      update_pull_coord_reference_value(pcrd, coord_ind, t);
  
+     get_pull_coord_distance(pull, coord_ind, pbc);
      /* Determine the deviation */
      dev = pcrd->value - pcrd->value_ref;
  
@@@ -1674,7 -1677,7 +1677,7 @@@ void apply_external_pull_coord_force(st
          calc_pull_coord_vector_force(pcrd);
  
          /* Add the forces for this coordinate to the total virial and force */
-         if (forceWithVirial->computeVirial_)
+         if (forceWithVirial->computeVirial_ && pull->comm.isMasterRank)
          {
              matrix virial = { { 0 } };
              add_virial_coord(virial, pcrd);
@@@ -1710,7 -1713,7 +1713,7 @@@ static void do_pull_pot_coord(struct pu
  }
  
  real pull_potential(struct pull_t *pull, t_mdatoms *md, t_pbc *pbc,
 -                    t_commrec *cr, double t, real lambda,
 +                    const t_commrec *cr, double t, real lambda,
                      rvec *x, gmx::ForceWithVirial *force, real *dvdlambda)
  {
      real V = 0;
  }
  
  void pull_constraint(struct pull_t *pull, t_mdatoms *md, t_pbc *pbc,
 -                     t_commrec *cr, double dt, double t,
 +                     const t_commrec *cr, double dt, double t,
                       rvec *x, rvec *xp, rvec *v, tensor vir)
  {
      assert(pull != NULL);
@@@ -1817,7 -1820,7 +1820,7 @@@ static void make_local_pull_group(gmx_g
      }
  }
  
 -void dd_make_local_pull_groups(t_commrec *cr, struct pull_t *pull, t_mdatoms *md)
 +void dd_make_local_pull_groups(const t_commrec *cr, struct pull_t *pull, t_mdatoms *md)
  {
      gmx_domdec_t   *dd;
      pull_comm_t    *comm;
          ga2la = nullptr;
      }
  
-     /* We always make the master node participate, such that it can do i/o
-      * and to simplify MC type extensions people might have.
+     /* We always make the master node participate, such that it can do i/o,
+      * add the virial and to simplify MC type extensions people might have.
       */
-     bMustParticipate = (comm->bParticipateAll || dd == nullptr || DDMASTER(dd));
+     bMustParticipate = (comm->bParticipateAll || comm->isMasterRank);
  
      for (g = 0; g < pull->ngroup; g++)
      {
          make_local_pull_group(ga2la, &pull->group[g],
                                0, md->homenr);
  
+         GMX_ASSERT(bMustParticipate || dd != nullptr, "Either all ranks (including this rank) participate, or we use DD and need to have access to dd here");
          /* We should participate if we have pull or pbc atoms */
          if (!bMustParticipate &&
              (pull->group[g].nat_loc > 0 ||
      pull->bSetPBCatoms = TRUE;
  }
  
 -static void init_pull_group_index(FILE *fplog, t_commrec *cr,
 +static void init_pull_group_index(FILE *fplog, const t_commrec *cr,
                                    int g, pull_group_work_t *pg,
                                    gmx_bool bConstraint, ivec pulldim_con,
                                    const gmx_mtop_t *mtop,
  struct pull_t *
  init_pull(FILE *fplog, const pull_params_t *pull_params, const t_inputrec *ir,
            int nfile, const t_filenm fnm[],
 -          const gmx_mtop_t *mtop, t_commrec *cr,
 +          const gmx_mtop_t *mtop, const t_commrec *cr,
            const gmx_output_env_t *oenv, real lambda,
            gmx_bool bOutFile,
            const ContinuationOptions &continuationOptions)
      comm = &pull->comm;
  
  #if GMX_MPI
-     /* Use a sub-communicator when we have more than 32 ranks */
+     /* Use a sub-communicator when we have more than 32 ranks, but not
+      * when we have an external pull potential, since then the external
+      * potential provider expects each rank to have the coordinate.
+      */
      comm->bParticipateAll = (cr == nullptr || !DOMAINDECOMP(cr) ||
                               cr->dd->nnodes <= 32 ||
+                              pull->numCoordinatesWithExternalPotential > 0 ||
                               getenv("GMX_PULL_PARTICIPATE_ALL") != nullptr);
      /* This sub-commicator is not used with comm->bParticipateAll,
       * so we can always initialize it to NULL.
       */
      comm->mpi_comm_com    = MPI_COMM_NULL;
      comm->nparticipate    = 0;
+     comm->isMasterRank    = (cr == nullptr || MASTER(cr));
  #else
      /* No MPI: 1 rank: all ranks pull */
      comm->bParticipateAll = TRUE;
+     comm->isMasterRank    = true;
  #endif
      comm->bParticipate    = comm->bParticipateAll;
      comm->setup_count     = 0;
index d302c1f425b37a62b3eb352e2802eb867e46dedb,1695e0b72d53729f4ece973596131895d9a122b8..a8319e4bf8c89905242648c97b6766d92a4ec039
@@@ -95,6 -95,7 +95,7 @@@
  #include "gromacs/mdlib/vcm.h"
  #include "gromacs/mdlib/vsite.h"
  #include "gromacs/mdtypes/awh-history.h"
+ #include "gromacs/mdtypes/awh-params.h"
  #include "gromacs/mdtypes/commrec.h"
  #include "gromacs/mdtypes/df_history.h"
  #include "gromacs/mdtypes/energyhistory.h"
@@@ -280,9 -281,7 +281,9 @@@ static void prepareRerunState(const t_t
  }
  
  /*! \libinternal
 -    \copydoc integrator_t (FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
 +    \copydoc integrator_t (FILE *fplog, t_commrec *cr,
 +                           const gmx_multisim_t *ms,
 +                           const gmx::MDLogger &mdlog,
                             int nfile, const t_filenm fnm[],
                             const gmx_output_env_t *oenv,
                             const MdrunOptions &mdrunOptions,
                             gmx_membed_t *membed,
                             gmx_walltime_accounting_t walltime_accounting)
   */
 -double gmx::do_md(FILE *fplog, t_commrec *cr, const gmx::MDLogger &mdlog,
 +double gmx::do_md(FILE *fplog, t_commrec *cr,
 +                  const gmx_multisim_t *ms,
 +                  const gmx::MDLogger &mdlog,
                    int nfile, const t_filenm fnm[],
                    const gmx_output_env_t *oenv,
                    const MdrunOptions &mdrunOptions,
                    gmx_membed_t *membed,
                    gmx_walltime_accounting_t walltime_accounting)
  {
--    gmx_mdoutf_t    outf = nullptr;
--    gmx_int64_t     step, step_rel;
--    double          elapsed_time;
--    double          t, t0, lam0[efptNR];
--    gmx_bool        bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
--    gmx_bool        bNS, bNStList, bSimAnn, bStopCM,
--                    bFirstStep, bInitStep, bLastStep = FALSE,
-                     bUsingEnsembleRestraints;
 -                    bBornRadii;
++    gmx_mdoutf_t      outf = nullptr;
++    gmx_int64_t       step, step_rel;
++    double            elapsed_time;
++    double            t, t0, lam0[efptNR];
++    gmx_bool          bGStatEveryStep, bGStat, bCalcVir, bCalcEnerStep, bCalcEner;
++    gmx_bool          bNS, bNStList, bSimAnn, bStopCM,
++                      bFirstStep, bInitStep, bLastStep = FALSE;
      gmx_bool          bDoDHDL = FALSE, bDoFEP = FALSE, bDoExpanded = FALSE;
      gmx_bool          do_ene, do_log, do_verbose, bRerunWarnNoV = TRUE,
                        bForceUpdate = FALSE, bCPT;
      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);
 +    SimulationSignaller nullSignaller(nullptr, nullptr, nullptr, false, false);
  
      if (!mdrunOptions.writeConfout)
      {
      }
  
      /* Set up interactive MD (IMD) */
 -    init_IMD(ir, cr, top_global, fplog, ir->nstcalcenergy, MASTER(cr) ? as_rvec_array(state_global->x.data()) : nullptr,
 +    init_IMD(ir, cr, ms, top_global, fplog, ir->nstcalcenergy,
 +             MASTER(cr) ? as_rvec_array(state_global->x.data()) : nullptr,
               nfile, fnm, oenv, mdrunOptions);
  
      if (DOMAINDECOMP(cr))
      /* Initialize AWH and restore state from history in checkpoint if needed. */
      if (ir->bDoAwh)
      {
 -        ir->awh = new gmx::Awh(fplog, *ir, cr, *ir->awhParams, opt2fn("-awh", nfile, fnm), ir->pull_work);
 +        ir->awh = new gmx::Awh(fplog, *ir, cr, ms, *ir->awhParams, opt2fn("-awh", nfile, fnm), ir->pull_work);
  
          if (startingFromCheckpoint)
          {
      const bool useReplicaExchange = (replExParams.exchangeInterval > 0);
      if (useReplicaExchange && MASTER(cr))
      {
 -        repl_ex = init_replica_exchange(fplog, cr->ms, top_global->natoms, ir,
 +        repl_ex = init_replica_exchange(fplog, ms, top_global->natoms, ir,
                                          replExParams);
      }
--
      /* PME tuning is only supported in the Verlet scheme, with PME for
       * Coulomb. It is not supported with only LJ PME, or for
       * reruns. */
          {
              /* Constrain the initial coordinates and velocities */
              do_constrain_first(fplog, constr, ir, mdatoms, state,
 -                               cr, nrnb, fr, top);
 +                               cr, ms, nrnb, fr, top);
          }
          if (vsite)
          {
      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) || (ms && fcd->orires.nr);
  
+     bool simulationsShareState = false;
+     int  nstSignalComm         = nstglobalcomm;
      {
-         // Replica exchange and ensemble restraints need all
+         // TODO This implementation of ensemble orientation restraints is nasty because
+         // a user can't just do multi-sim with single-sim orientation restraints.
 -        bool usingEnsembleRestraints = (fcd->disres.nsystems > 1) || (cr->ms && fcd->orires.nr);
 -        bool awhUsesMultiSim         = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && cr->ms);
++        bool usingEnsembleRestraints = (fcd->disres.nsystems > 1) || (ms && fcd->orires.nr);
++        bool awhUsesMultiSim         = (ir->bDoAwh && ir->awhParams->shareBiasMultisim && ms);
+         // Replica exchange, ensemble restraints and AWH 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    = !useReplicaExchange && !bUsingEnsembleRestraints;
-         bool stopConditionIsLocal = !useReplicaExchange && !bUsingEnsembleRestraints;
+         // TODO: Make algorithm initializers set these flags.
+         simulationsShareState     = useReplicaExchange || usingEnsembleRestraints || awhUsesMultiSim;
          bool resetCountersIsLocal = true;
-         signals[eglsCHKPT]         = SimulationSignal(checkpointIsLocal);
-         signals[eglsSTOPCOND]      = SimulationSignal(stopConditionIsLocal);
+         signals[eglsCHKPT]         = SimulationSignal(!simulationsShareState);
+         signals[eglsSTOPCOND]      = SimulationSignal(!simulationsShareState);
          signals[eglsRESETCOUNTERS] = SimulationSignal(resetCountersIsLocal);
+         if (simulationsShareState)
+         {
+             // Inter-simulation signal communication does not need to happen
+             // often, so we use a minimum of 200 steps to reduce overhead.
+             const int c_minimumInterSimulationSignallingInterval = 200;
+             nstSignalComm = ((c_minimumInterSimulationSignallingInterval + nstglobalcomm - 1)/nstglobalcomm)*nstglobalcomm;
+         }
      }
  
      DdOpenBalanceRegionBeforeForceComputation ddOpenBalanceRegion   = (DOMAINDECOMP(cr) ? DdOpenBalanceRegionBeforeForceComputation::yes : DdOpenBalanceRegionBeforeForceComputation::no);
      step_rel = 0;
  
      // TODO extract this to new multi-simulation module
 -    if (MASTER(cr) && MULTISIM(cr) && !useReplicaExchange)
 +    if (MASTER(cr) && isMultiSim(ms) && !useReplicaExchange)
      {
 -        if (!multisim_int_all_are_equal(cr->ms, ir->nsteps))
 +        if (!multisim_int_all_are_equal(ms, ir->nsteps))
          {
              GMX_LOG(mdlog.warning).appendText(
                      "Note: The number of steps is not consistent across multi simulations,\n"
                      "but we are proceeding anyway!");
          }
 -        if (!multisim_int_all_are_equal(cr->ms, ir->init_step))
 +        if (!multisim_int_all_are_equal(ms, ir->init_step))
          {
              GMX_LOG(mdlog.warning).appendText(
                      "Note: The initial step is not consistent across multi simulations,\n"
              bLastStep = TRUE;
          }
  
 -        /* Determine whether or not to update the Born radii if doing GB */
 -        bBornRadii = bFirstStep;
 -        if (ir->implicit_solvent && (step % ir->nstgbradii == 0))
 -        {
 -            bBornRadii = TRUE;
 -        }
 -
          /* do_log triggers energy and virial calculation. Because this leads
           * to different code paths, forces can be different. Thus for exact
           * continuation we should avoid extra log output.
          if (shellfc)
          {
              /* Now is the time to relax the shells */
 -            relax_shell_flexcon(fplog, cr, mdrunOptions.verbose, step,
 +            relax_shell_flexcon(fplog, cr, ms, mdrunOptions.verbose, step,
                                  ir, bNS, force_flags, top,
                                  constr, enerd, fcd,
                                  state, &f, force_vir, mdatoms,
                                  nrnb, wcycle, graph, groups,
 -                                shellfc, fr, bBornRadii, t, mu_tot,
 +                                shellfc, fr, t, mu_tot,
                                  vsite,
                                  ddOpenBalanceRegion, ddCloseBalanceRegion);
          }
               * This is parallellized as well, and does communication too.
               * Check comments in sim_util.c
               */
 -            do_force(fplog, cr, ir, step, nrnb, wcycle, top, groups,
 +            do_force(fplog, cr, ms, ir, step, nrnb, wcycle, top, groups,
                       state->box, state->x, &state->hist,
                       f, force_vir, mdatoms, enerd, fcd,
                       state->lambda, graph,
 -                     fr, vsite, mu_tot, t, ed, bBornRadii,
 +                     fr, vsite, mu_tot, t, ed,
                       (bNS ? GMX_FORCE_NS : 0) | force_flags,
                       ddOpenBalanceRegion, ddCloseBalanceRegion);
          }
                  update_constraints(fplog, step, nullptr, ir, mdatoms,
                                     state, fr->bMolPBC, graph, f,
                                     &top->idef, shake_vir,
 -                                   cr, nrnb, wcycle, upd, constr,
 +                                   cr, ms, nrnb, wcycle, upd, constr,
                                     TRUE, bCalcVir);
                  wallcycle_start(wcycle, ewcUPDATE);
              }
                   * to allow for exact continuation, when possible.
                   */
                  signals[eglsSTOPCOND].sig = 1;
-                 nsteps_stop               = std::max(ir->nstlist, 2*nstglobalcomm);
+                 nsteps_stop               = std::max(ir->nstlist, 2*nstSignalComm);
              }
              else if (gmx_get_stop_condition() == gmx_stop_cond_next)
              {
                   * This breaks exact continuation.
                   */
                  signals[eglsSTOPCOND].sig = -1;
-                 nsteps_stop               = nstglobalcomm + 1;
+                 nsteps_stop               = nstSignalComm + 1;
              }
              if (fplog)
              {
                  update_constraints(fplog, step, nullptr, ir, mdatoms,
                                     state, fr->bMolPBC, graph, f,
                                     &top->idef, tmp_vir,
 -                                   cr, nrnb, wcycle, upd, constr,
 +                                   cr, ms, nrnb, wcycle, upd, constr,
                                     TRUE, bCalcVir);
              }
          }
              update_constraints(fplog, step, &dvdl_constr, ir, mdatoms, state,
                                 fr->bMolPBC, graph, f,
                                 &top->idef, shake_vir,
 -                               cr, nrnb, wcycle, upd, constr,
 +                               cr, ms, nrnb, wcycle, upd, constr,
                                 FALSE, bCalcVir);
  
              if (ir->eI == eiVVAK)
                  update_constraints(fplog, step, nullptr, ir, mdatoms,
                                     state, fr->bMolPBC, graph, f,
                                     &top->idef, tmp_vir,
 -                                   cr, nrnb, wcycle, upd, nullptr,
 +                                   cr, ms, nrnb, wcycle, upd, nullptr,
                                     FALSE, bCalcVir);
              }
              if (EI_VV(ir->eI))
          {
              // Organize to do inter-simulation signalling on steps if
              // and when algorithms require it.
-             bool doInterSimSignal = (!bFirstStep && bDoReplEx) || bUsingEnsembleRestraints;
+             bool doInterSimSignal = (simulationsShareState && do_per_step(step, nstSignalComm));
  
              if (bGStat || (!EI_VV(ir->eI) && do_per_step(step+1, nstglobalcomm)) || doInterSimSignal)
              {
                  // situation where e.g. checkpointing can't be
                  // signalled.
                  bool                doIntraSimSignal = true;
 -                SimulationSignaller signaller(&signals, cr, doInterSimSignal, doIntraSimSignal);
 +                SimulationSignaller signaller(&signals, cr, ms, 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,
              state->fep_state = lamnew;
          }
          /* Print the remaining wall clock time for the run */
 -        if (MULTIMASTER(cr) &&
 +        if (isMasterSimMasterRank(ms, cr) &&
              (do_verbose || gmx_got_usr_signal()) &&
              !bPMETunePrinting)
          {
          bExchanged = FALSE;
          if (bDoReplEx)
          {
 -            bExchanged = replica_exchange(fplog, cr, repl_ex,
 +            bExchanged = replica_exchange(fplog, cr, ms, repl_ex,
                                            state_global, enerd,
                                            state, step, t);
          }
                         eprAVER, mdebin, fcd, groups, &(ir->opts), ir->awh);
          }
      }
 -    // TODO Enable the next line when merging to master branch
 -    // done_ebin(mdebin->ebin);
 +    done_ebin(mdebin->ebin);
      done_mdoutf(outf);
  
      if (bPMETune)