Merge branch release-2018 into release-2019
[alexxy/gromacs.git] / src / gromacs / fileio / checkpoint.cpp
index ee267914eec05741972e2d8b2a7c081d360a7dea..2dcf992ffc2cc532fe2f1805b36b8a4b6675ce0e 100644 (file)
 #include <sys/locking.h>
 #endif
 
+#include <array>
+#include <memory>
+
 #include "buildinfo.h"
+#include "gromacs/compat/make_unique.h"
 #include "gromacs/fileio/filetypes.h"
 #include "gromacs/fileio/gmxfio.h"
 #include "gromacs/fileio/gmxfio-xdr.h"
@@ -70,6 +74,7 @@
 #include "gromacs/mdtypes/inputrec.h"
 #include "gromacs/mdtypes/md_enums.h"
 #include "gromacs/mdtypes/observableshistory.h"
+#include "gromacs/mdtypes/pullhistory.h"
 #include "gromacs/mdtypes/state.h"
 #include "gromacs/mdtypes/swaphistory.h"
 #include "gromacs/trajectory/trajectoryframe.h"
@@ -85,7 +90,7 @@
 #include "gromacs/utility/sysinfo.h"
 #include "gromacs/utility/txtdump.h"
 
-#ifdef GMX_FAHCORE
+#if GMX_FAHCORE
 #include "corewrap.h"
 #endif
 
 #define CPT_MAGIC2 171819
 #define CPTSTRLEN 1024
 
-/* cpt_version should normally only be changed
- * when the header or footer format changes.
+/*! \brief Enum of values that describe the contents of a cpt file
+ * whose format matches a version number
+ *
+ * The enum helps the code be more self-documenting and ensure merges
+ * do not silently resolve when two patches make the same bump. When
+ * adding new functionality, add a new element just above cptv_Count
+ * in this enumeration, and write code below that does the right thing
+ * according to the value of file_version.
+ */
+enum cptv {
+    cptv_Unknown = 17,                                       /**< Version before numbering scheme */
+    cptv_RemoveBuildMachineInformation,                      /**< remove functionality that makes mdrun builds non-reproducible */
+    cptv_ComPrevStepAsPullGroupReference,                    /**< Allow using COM of previous step as pull group PBC reference */
+    cptv_PullAverage,                                        /**< Added possibility to output average pull force and position */
+    cptv_Count                                               /**< the total number of cptv versions */
+};
+
+/*! \brief Version number of the file format written to checkpoint
+ * files by this version of the code.
+ *
+ * cpt_version should normally only be changed, via adding a new field
+ * to cptv enumeration, when the header or footer format changes.
+ *
  * The state data format itself is backward and forward compatible.
  * But old code can not read a new entry that is present in the file
  * (but can read a new format when new entries are not present).
- */
-static const int cpt_version = 17;
+ *
+ * The cpt_version increases whenever the file format in the main
+ * development branch changes, due to an extension of the cptv enum above.
+ * Backward compatibility for reading old run input files is maintained
+ * by checking this version number against that of the file and then using
+ * the correct code path. */
+static const int cpt_version = cptv_Count - 1;
 
 
 const char *est_names[estNR] =
@@ -110,7 +141,7 @@ const char *est_names[estNR] =
     "x", "v", "sdx-unsupported", "CGp", "LD-rng-unsupported", "LD-rng-i-unsupported",
     "disre_initf", "disre_rm3tav",
     "orire_initf", "orire_Dtav",
-    "svir_prev", "nosehoover-vxi", "v_eta", "vol0", "nhpres_xi", "nhpres_vxi", "fvir_prev", "fep_state", "MC-rng-unsupported", "MC-rng-i-unsupported"
+    "svir_prev", "nosehoover-vxi", "v_eta", "vol0", "nhpres_xi", "nhpres_vxi", "fvir_prev", "fep_state", "MC-rng-unsupported", "MC-rng-i-unsupported",
     "barostat-integral"
 };
 
@@ -118,7 +149,7 @@ enum {
     eeksEKIN_N, eeksEKINH, eeksDEKINDL, eeksMVCOS, eeksEKINF, eeksEKINO, eeksEKINSCALEF, eeksEKINSCALEH, eeksVSCALE, eeksEKINTOTAL, eeksNR
 };
 
-const char *eeks_names[eeksNR] =
+static const char *eeks_names[eeksNR] =
 {
     "Ekin_n", "Ekinh", "dEkindlambda", "mv_cos",
     "Ekinf", "Ekinh_old", "EkinScaleF_NHC", "EkinScaleH_NHC", "Vscale_NHC", "Ekin_Total"
@@ -135,7 +166,22 @@ enum {
     eenhNR
 };
 
-const char *eenh_names[eenhNR] =
+enum {
+    epullhPULL_NUMCOORDINATES, epullhPULL_NUMGROUPS, epullhPULL_NUMVALUESINXSUM,
+    epullhPULL_NUMVALUESINFSUM,  epullhNR
+};
+
+enum {
+    epullcoordh_VALUE_REF_SUM, epullcoordh_VALUE_SUM, epullcoordh_DR01_SUM,
+    epullcoordh_DR23_SUM, epullcoordh_DR45_SUM, epullcoordh_FSCAL_SUM,
+    epullcoordh_DYNAX_SUM, epullcoordh_NR
+};
+
+enum {
+    epullgrouph_X_SUM, epullgrouph_NR
+};
+
+static const char *eenh_names[eenhNR] =
 {
     "energy_n", "energy_aver", "energy_sum", "energy_nsum",
     "energy_sum_sim", "energy_nsum_sim",
@@ -146,13 +192,19 @@ const char *eenh_names[eenhNR] =
     "energy_delta_h_start_lambda"
 };
 
+static const char *ePullhNames[epullhNR] =
+{
+    "pullhistory_numcoordinates", "pullhistory_numgroups", "pullhistory_numvaluesinxsum",
+    "pullhistory_numvaluesinfsum"
+};
+
 /* free energy history variables -- need to be preserved over checkpoint */
 enum {
     edfhBEQUIL, edfhNATLAMBDA, edfhWLHISTO, edfhWLDELTA, edfhSUMWEIGHTS, edfhSUMDG, edfhSUMMINVAR, edfhSUMVAR,
     edfhACCUMP, edfhACCUMM, edfhACCUMP2, edfhACCUMM2, edfhTIJ, edfhTIJEMP, edfhNR
 };
 /* free energy history variable names  */
-const char *edfh_names[edfhNR] =
+static const char *edfh_names[edfhNR] =
 {
     "bEquilibrated", "N_at_state", "Wang-Landau Histogram", "Wang-Landau Delta", "Weights", "Free Energies", "minvar", "variance",
     "accumulated_plus", "accumulated_minus", "accumulated_plus_2",  "accumulated_minus_2", "Tij", "Tij_empirical"
@@ -172,7 +224,7 @@ enum {
     eawhhNR
 };
 
-const char *eawhh_names[eawhhNR] =
+static const char *eawhh_names[eawhhNR] =
 {
     "awh_in_initial",
     "awh_equilibrateHistogram",
@@ -185,6 +237,17 @@ const char *eawhh_names[eawhhNR] =
     "awh_forceCorrelationGrid"
 };
 
+enum {
+    epullsPREVSTEPCOM,
+    epullsNR
+};
+
+static const char *epull_prev_step_com_names[epullsNR] =
+{
+    "Pull groups prev step COM"
+};
+
+
 //! Higher level vector element type, only used for formatting checkpoint dumps
 enum class CptElementType
 {
@@ -201,7 +264,9 @@ enum class StatePart
     kineticEnergy,      //!< Kinetic energy, needed for T/P-coupling state
     energyHistory,      //!< Energy observable statistics
     freeEnergyHistory,  //!< Free-energy state and observable statistics
-    accWeightHistogram  //!< Accelerated weight histogram method state
+    accWeightHistogram, //!< Accelerated weight histogram method state
+    pullState,          //!< COM of previous step.
+    pullHistory         //!< Pull history statistics (sums since last written output)
 };
 
 //! \brief Return the name of a checkpoint entry based on part and part entry
@@ -214,6 +279,8 @@ static const char *entryName(StatePart part, int ecpt)
         case StatePart::energyHistory:      return eenh_names[ecpt];
         case StatePart::freeEnergyHistory:  return edfh_names[ecpt];
         case StatePart::accWeightHistogram: return eawhh_names[ecpt];
+        case StatePart::pullState:          return epull_prev_step_com_names[ecpt];
+        case StatePart::pullHistory:        return ePullhNames[ecpt];
     }
 
     return nullptr;
@@ -224,25 +291,21 @@ static void cp_warning(FILE *fp)
     fprintf(fp, "\nWARNING: Checkpoint file is corrupted or truncated\n\n");
 }
 
-static void cp_error()
+[[ noreturn ]] 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);
     }
 }
 
@@ -294,7 +357,19 @@ static void do_cpt_int_err(XDR *xd, const char *desc, int *i, FILE *list)
     }
 }
 
-static void do_cpt_step_err(XDR *xd, const char *desc, gmx_int64_t *i, FILE *list)
+static void do_cpt_bool_err(XDR *xd, const char *desc, bool *b, FILE *list)
+{
+    int i   = static_cast<int>(*b);
+
+    if (do_cpt_int(xd, desc, &i, list) < 0)
+    {
+        cp_error();
+    }
+
+    *b = (i != 0);
+}
+
+static void do_cpt_step_err(XDR *xd, const char *desc, int64_t *i, FILE *list)
 {
     char   buf[STEPSTRSIZE];
 
@@ -357,21 +432,18 @@ struct xdr_type
 template <>
 struct xdr_type<int>
 {
-    // cppcheck-suppress unusedStructMember
     static const int value = xdr_datatype_int;
 };
 
 template <>
 struct xdr_type<float>
 {
-    // cppcheck-suppress unusedStructMember
     static const int value = xdr_datatype_float;
 };
 
 template <>
 struct xdr_type<double>
 {
-    // cppcheck-suppress unusedStructMember
     static const int value = xdr_datatype_double;
 };
 
@@ -382,13 +454,10 @@ static inline unsigned int sizeOfXdrType(int xdrType)
     {
         case xdr_datatype_int:
             return sizeof(int);
-            break;
         case xdr_datatype_float:
             return sizeof(float);
-            break;
         case xdr_datatype_double:
             return sizeof(double);
-            break;
         default: GMX_RELEASE_ASSERT(false, "XDR data type not implemented");
     }
 
@@ -402,13 +471,10 @@ static inline xdrproc_t xdrProc(int xdrType)
     {
         case xdr_datatype_int:
             return reinterpret_cast<xdrproc_t>(xdr_int);
-            break;
         case xdr_datatype_float:
             return reinterpret_cast<xdrproc_t>(xdr_float);
-            break;
         case xdr_datatype_double:
             return reinterpret_cast<xdrproc_t>(xdr_double);
-            break;
         default: GMX_RELEASE_ASSERT(false, "XDR data type not implemented");
     }
 
@@ -442,14 +508,12 @@ static bool_t listXdrVector(XDR *xd, StatePart part, int ecpt, int nf, int xdrTy
 #if !GMX_DOUBLE
                 if (cptElementType == CptElementType::real3)
                 {
-                    // cppcheck-suppress invalidPointerCast
                     pr_rvecs(list, 0, entryName(part, ecpt), reinterpret_cast<const rvec *>(data.data()), nf/3);
                 }
                 else
 #endif
                 {
                     /* Note: With double precision code dumping a single precision rvec will produce float iso rvec print, but that's a minor annoyance */
-                    // cppcheck-suppress invalidPointerCast
                     pr_fvec(list, 0, entryName(part, ecpt), reinterpret_cast<const float *>(data.data()), nf, TRUE);
                 }
                 break;
@@ -457,14 +521,12 @@ static bool_t listXdrVector(XDR *xd, StatePart part, int ecpt, int nf, int xdrTy
 #if GMX_DOUBLE
                 if (cptElementType == CptElementType::real3)
                 {
-                    // cppcheck-suppress invalidPointerCast
                     pr_rvecs(list, 0, entryName(part, ecpt), reinterpret_cast<const rvec *>(data.data()), nf/3);
                 }
                 else
 #endif
                 {
                     /* Note: With single precision code dumping a double precision rvec will produce float iso rvec print, but that's a minor annoyance */
-                    // cppcheck-suppress invalidPointerCast
                     pr_dvec(list, 0, entryName(part, ecpt), reinterpret_cast<const double *>(data.data()), nf, TRUE);
                 }
                 break;
@@ -478,7 +540,6 @@ static bool_t listXdrVector(XDR *xd, StatePart part, int ecpt, int nf, int xdrTy
 //! \brief Convert a double array, typed char*, to float
 gmx_unused static void convertArrayRealPrecision(const char *c, float *v, int n)
 {
-    // cppcheck-suppress invalidPointerCast
     const double *d = reinterpret_cast<const double *>(c);
     for (int i = 0; i < n; i++)
     {
@@ -489,7 +550,6 @@ gmx_unused static void convertArrayRealPrecision(const char *c, float *v, int n)
 //! \brief Convert a float array, typed char*, to double
 static void convertArrayRealPrecision(const char *c, double *v, int n)
 {
-    // cppcheck-suppress invalidPointerCast
     const float *f = reinterpret_cast<const float *>(c);
     for (int i = 0; i < n; i++)
     {
@@ -545,7 +605,6 @@ static int doVectorLow(XDR *xd, StatePart part, int ecpt, int sflags,
             if (v != nullptr)
             {
                 GMX_RELEASE_ASSERT(nptr != nullptr, "With nval<0 we should have nptr!=NULL");
-                // cppcheck-suppress nullPointer
                 numElemInTheFile = *nptr;
             }
             else
@@ -598,7 +657,6 @@ static int doVectorLow(XDR *xd, StatePart part, int ecpt, int sflags,
             {
                 gmx_fatal(FARGS, "Type %s: incompatible checkpoint formats or corrupted checkpoint file.", buf);
             }
-            fprintf(stderr, "Precision %s\n", buf);
         }
 
         T *vp;
@@ -614,7 +672,7 @@ static int doVectorLow(XDR *xd, StatePart part, int ecpt, int sflags,
         {
             /* This conditional ensures that we don't resize on write.
              * In particular in the state where this code was written
-             * PaddedRVecVector has a size of numElemInThefile and we
+             * vector has a size of numElemInThefile and we
              * don't want to lose that padding here.
              */
             if (vector->size() < static_cast<unsigned int>(numElemInTheFile))
@@ -661,12 +719,12 @@ static int doVectorLow(XDR *xd, StatePart part, int ecpt, int sflags,
     return 0;
 }
 
-//! \brief Read/Write a std::vector.
+//! \brief Read/Write a std::vector, on read checks the number of elements matches \p numElements, if specified.
 template <typename T>
 static int doVector(XDR *xd, StatePart part, int ecpt, int sflags,
-                    std::vector<T> *vector, FILE *list)
+                    std::vector<T> *vector, FILE *list, int numElements = -1)
 {
-    return doVectorLow<T>(xd, part, ecpt, sflags, -1, nullptr, nullptr, vector, list, CptElementType::real);
+    return doVectorLow<T>(xd, part, ecpt, sflags, numElements, nullptr, nullptr, vector, list, CptElementType::real);
 }
 
 //! \brief Read/Write an ArrayRef<real>.
@@ -677,31 +735,32 @@ static int doRealArrayRef(XDR *xd, StatePart part, int ecpt, int sflags,
     return doVectorLow < real, std::allocator < real>>(xd, part, ecpt, sflags, vector.size(), nullptr, &v_real, nullptr, list, CptElementType::real);
 }
 
-//! \brief Read/Write a vector whose value_type is RVec and whose allocator is \c AllocatorType.
-template <typename AllocatorType>
+//! Convert from view of RVec to view of real.
+static gmx::ArrayRef<real>
+realArrayRefFromRVecArrayRef(gmx::ArrayRef<gmx::RVec> ofRVecs)
+{
+    return gmx::arrayRefFromArray<real>(reinterpret_cast<real *>(ofRVecs.data()), ofRVecs.size() * DIM);
+}
+
+//! \brief Read/Write a PaddedVector whose value_type is RVec.
+template <typename PaddedVectorOfRVecType>
 static int doRvecVector(XDR *xd, StatePart part, int ecpt, int sflags,
-                        std::vector<gmx::RVec, AllocatorType> *v, int numAtoms, FILE *list)
+                        PaddedVectorOfRVecType *v, int numAtoms, FILE *list)
 {
     const int numReals = numAtoms*DIM;
 
     if (list == nullptr)
     {
         GMX_RELEASE_ASSERT(sflags & (1 << ecpt), "When not listing, the flag for the entry should be set when requesting i/o");
-        GMX_RELEASE_ASSERT(v->size() >= static_cast<size_t>(numAtoms), "v should have sufficient size for numAtoms");
-
-        real *v_real = v->data()->as_vec();
+        GMX_RELEASE_ASSERT(v->size() == numAtoms, "v should have sufficient size for numAtoms");
 
-        // PaddedRVecVector is padded beyond numAtoms, we should only write
-        // numAtoms RVecs
-        gmx::ArrayRef<real> ref(v_real, v_real + numReals);
-
-        return doRealArrayRef(xd, part, ecpt, sflags, ref, list);
+        return doRealArrayRef(xd, part, ecpt, sflags, realArrayRefFromRVecArrayRef(makeArrayRef(*v)), list);
     }
     else
     {
         // Use the rebind facility to change the value_type of the
         // allocator from RVec to real.
-        using realAllocator = typename AllocatorType::template rebind<real>::other;
+        using realAllocator = typename std::allocator_traits<typename PaddedVectorOfRVecType::allocator_type>::template rebind_alloc<real>;
         return doVectorLow<real, realAllocator>(xd, part, ecpt, sflags, numReals, nullptr, nullptr, nullptr, list, CptElementType::real);
     }
 }
@@ -744,6 +803,15 @@ static int do_cpte_int(XDR *xd, StatePart part, int ecpt, int sflags,
     return do_cpte_ints(xd, part, ecpt, sflags, 1, &i, list);
 }
 
+static int do_cpte_bool(XDR *xd, StatePart part, int ecpt, int sflags,
+                        bool *b, FILE *list)
+{
+    int i   = static_cast<int>(*b);
+    int ret = do_cpte_int(xd, part, ecpt, sflags, &i, list);
+    *b = (i != 0);
+    return ret;
+}
+
 static int do_cpte_doubles(XDR *xd, StatePart part, int ecpt, int sflags,
                            int n, double **v, FILE *list)
 {
@@ -876,23 +944,55 @@ static int do_cpte_matrices(XDR *xd, StatePart part, int ecpt, int sflags,
     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.
+//
+// Several fields were once written to checkpoint file headers, but
+// have been removed. So that old files can continue to be read,
+// the names of such fields contain the string "_UNUSED" so that it
+// is clear they should not be used.
+struct CheckpointHeaderContents
+{
+    int         file_version;
+    char        version[CPTSTRLEN];
+    char        btime_UNUSED[CPTSTRLEN];
+    char        buser_UNUSED[CPTSTRLEN];
+    char        bhost_UNUSED[CPTSTRLEN];
+    int         double_prec;
+    char        fprog[CPTSTRLEN];
+    char        ftime[CPTSTRLEN];
+    int         eIntegrator;
+    int         simulation_part;
+    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         flagsPullHistory;
+    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)
     {
@@ -913,138 +1013,149 @@ static void do_cpt_header(XDR *xd, gmx_bool bRead, int *file_version,
                   "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);
+    // The following fields are no longer ever written with meaningful
+    // content, but because they precede the file version, there is no
+    // good way for new code to read the old and new formats, nor a
+    // good way for old code to avoid giving an error while reading a
+    // new format. So we read and write a field that no longer has a
+    // purpose.
+    do_cpt_string_err(xd, "GROMACS build time UNUSED", contents->btime_UNUSED, list);
+    do_cpt_string_err(xd, "GROMACS build user UNUSED", contents->buser_UNUSED, list);
+    do_cpt_string_err(xd, "GROMACS build host UNUSED", contents->bhost_UNUSED, 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<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;
+    }
+
+    if (contents->file_version >= 18)
+    {
+        do_cpt_int_err(xd, "pull history flags", &contents->flagsPullHistory, list);
+    }
+    else
+    {
+        contents->flagsPullHistory = 0;
     }
 }
 
@@ -1116,13 +1227,13 @@ static int do_cpt_state(XDR *xd,
                 case estDISRE_RM3TAV: ret         = do_cpte_n_reals(xd, part, i, sflags, &state->hist.ndisrepairs, &state->hist.disre_rm3tav, list); break;
                 case estORIRE_INITF:  ret         = do_cpte_real (xd, part, i, sflags, &state->hist.orire_initf, list); break;
                 case estORIRE_DTAV:   ret         = do_cpte_n_reals(xd, part, i, sflags, &state->hist.norire_Dtav, &state->hist.orire_Dtav, list); break;
+                case estPREVSTEPCOM:  ret         = doVector<double>(xd, part, i, sflags, &state->com_prev_step, list); break;
                 default:
                     gmx_fatal(FARGS, "Unknown state entry %d\n"
                               "You are reading a checkpoint file written by different code, which is not supported", i);
             }
         }
     }
-
     return ret;
 }
 
@@ -1354,7 +1465,7 @@ static int do_cpt_enerhist(XDR *xd, gmx_bool bRead,
                     {
                         if (deltaH == nullptr)
                         {
-                            enerhist->deltaHForeignLambdas.reset(new delta_h_history_t);
+                            enerhist->deltaHForeignLambdas = gmx::compat::make_unique<delta_h_history_t>();
                             deltaH = enerhist->deltaHForeignLambdas.get();
                         }
                         deltaH->dh.resize(numDeltaH);
@@ -1401,6 +1512,140 @@ static int do_cpt_enerhist(XDR *xd, gmx_bool bRead,
     return ret;
 }
 
+static int doCptPullCoordHist(XDR *xd, PullCoordinateHistory *pullCoordHist,
+                              const StatePart part, FILE *list)
+{
+    int ret   = 0;
+    int flags = 0;
+
+    flags |= ((1<<epullcoordh_VALUE_REF_SUM) | (1<<epullcoordh_VALUE_SUM) | (1<<epullcoordh_DR01_SUM) |
+              (1<<epullcoordh_DR23_SUM) | (1<<epullcoordh_DR45_SUM) | (1<<epullcoordh_FSCAL_SUM));
+
+    for (int i = 0; i < epullcoordh_NR && ret == 0; i++)
+    {
+        switch (i)
+        {
+            case epullcoordh_VALUE_REF_SUM: ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->valueRef), list); break;
+            case epullcoordh_VALUE_SUM:     ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->value), list); break;
+            case epullcoordh_DR01_SUM:
+                for (int j = 0; j < DIM && ret == 0; j++)
+                {
+                    ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dr01[j]), list);
+                }
+                break;
+            case epullcoordh_DR23_SUM:
+                for (int j = 0; j < DIM && ret == 0; j++)
+                {
+                    ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dr23[j]), list);
+                }
+                break;
+            case epullcoordh_DR45_SUM:
+                for (int j = 0; j < DIM && ret == 0; j++)
+                {
+                    ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dr45[j]), list);
+                }
+                break;
+            case epullcoordh_FSCAL_SUM:     ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->scalarForce), list); break;
+            case epullcoordh_DYNAX_SUM:
+                for (int j = 0; j < DIM && ret == 0; j++)
+                {
+                    ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dynaX[j]), list);
+                }
+                break;
+        }
+    }
+
+    return ret;
+}
+
+static int doCptPullGroupHist(XDR *xd, PullGroupHistory *pullGroupHist,
+                              const StatePart part, FILE *list)
+{
+    int ret   = 0;
+    int flags = 0;
+
+    flags |= ((1<<epullgrouph_X_SUM));
+
+    for (int i = 0; i < epullgrouph_NR; i++)
+    {
+        switch (i)
+        {
+            case epullgrouph_X_SUM:
+                for (int j = 0; j < DIM && ret == 0; j++)
+                {
+                    ret = do_cpte_double(xd, part, i, flags, &(pullGroupHist->x[j]), list);
+                }
+                break;
+        }
+    }
+
+    return ret;
+}
+
+
+static int doCptPullHist(XDR *xd, gmx_bool bRead,
+                         int fflags, PullHistory *pullHist,
+                         const StatePart part,
+                         FILE *list)
+{
+    int ret                       = 0;
+    int pullHistoryNumCoordinates = 0;
+    int pullHistoryNumGroups      = 0;
+
+    /* Retain the number of terms in the sum and the number of coordinates (used for writing
+     * average pull forces and coordinates) in the pull history, in temporary variables,
+     * in case they cannot be read from the checkpoint, in order to have backward compatibility */
+    if (bRead)
+    {
+        pullHist->numValuesInXSum = 0;
+        pullHist->numValuesInFSum = 0;
+    }
+    else if (pullHist != nullptr)
+    {
+        pullHistoryNumCoordinates = pullHist->pullCoordinateSums.size();
+        pullHistoryNumGroups      = pullHist->pullGroupSums.size();
+    }
+    else
+    {
+        GMX_RELEASE_ASSERT(fflags == 0, "Without pull history, all flags should be off");
+    }
+
+    for (int i = 0; (i < epullhNR && ret == 0); i++)
+    {
+        if (fflags & (1<<i))
+        {
+            switch (i)
+            {
+                case epullhPULL_NUMCOORDINATES:         ret = do_cpte_int(xd, part, i, fflags, &pullHistoryNumCoordinates, list); break;
+                case epullhPULL_NUMGROUPS:              do_cpt_int_err(xd, eenh_names[i], &pullHistoryNumGroups, list); break;
+                case epullhPULL_NUMVALUESINXSUM:        do_cpt_int_err(xd, eenh_names[i], &pullHist->numValuesInXSum, list); break;
+                case epullhPULL_NUMVALUESINFSUM:        do_cpt_int_err(xd, eenh_names[i], &pullHist->numValuesInFSum, list); break;
+                default:
+                    gmx_fatal(FARGS, "Unknown pull history entry %d\n"
+                              "You are probably reading a new checkpoint file with old code", i);
+            }
+        }
+    }
+    if (bRead)
+    {
+        pullHist->pullCoordinateSums.resize(pullHistoryNumCoordinates);
+        pullHist->pullGroupSums.resize(pullHistoryNumGroups);
+    }
+    if (pullHist->numValuesInXSum > 0 || pullHist->numValuesInFSum > 0)
+    {
+        for (size_t i = 0; i < pullHist->pullCoordinateSums.size() && ret == 0; i++)
+        {
+            ret = doCptPullCoordHist(xd, &(pullHist->pullCoordinateSums[i]), part, list);
+        }
+        for (size_t i = 0; i < pullHist->pullGroupSums.size() && ret == 0; i++)
+        {
+            ret = doCptPullGroupHist(xd, &(pullHist->pullGroupSums[i]), part, list);
+        }
+    }
+
+    return ret;
+}
+
 static int do_cpt_df_hist(XDR *xd, int fflags, int nlambda, df_history_t **dfhistPtr, FILE *list)
 {
     int ret = 0;
@@ -1425,7 +1670,7 @@ static int do_cpt_df_hist(XDR *xd, int fflags, int nlambda, df_history_t **dfhis
         {
             switch (i)
             {
-                case edfhBEQUIL:       ret = do_cpte_int(xd, part, i, fflags, &dfhist->bEquil, list); break;
+                case edfhBEQUIL:       ret = do_cpte_bool(xd, part, i, fflags, &dfhist->bEquil, list); break;
                 case edfhNATLAMBDA:    ret = do_cpte_ints(xd, part, i, fflags, nlambda, &dfhist->n_at_lam, list); break;
                 case edfhWLHISTO:      ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->wl_histo, list); break;
                 case edfhWLDELTA:      ret = do_cpte_real(xd, part, i, fflags, &dfhist->wl_delta, list); break;
@@ -1559,9 +1804,9 @@ static int do_cpt_awh_bias(XDR *xd, gmx_bool bRead,
             switch (i)
             {
                 case eawhhIN_INITIAL:
-                    do_cpt_int_err(xd, eawhh_names[i], &state->in_initial, list); break;
+                    do_cpt_bool_err(xd, eawhh_names[i], &state->in_initial, list); break;
                 case eawhhEQUILIBRATEHISTOGRAM:
-                    do_cpt_int_err(xd, eawhh_names[i], &state->equilibrateHistogram, list); break;
+                    do_cpt_bool_err(xd, eawhh_names[i], &state->equilibrateHistogram, list); break;
                 case eawhhHISTSIZE:
                     do_cpt_double_err(xd, eawhh_names[i], &state->histSize, list); break;
                 case eawhhNPOINTS:
@@ -1633,7 +1878,7 @@ static int do_cpt_awh(XDR *xd, gmx_bool bRead,
         {
             GMX_RELEASE_ASSERT(bRead, "do_cpt_awh should not be called for writing without an AwhHistory");
 
-            awhHistoryLocal = std::shared_ptr<gmx::AwhHistory>(new gmx::AwhHistory());
+            awhHistoryLocal = std::make_shared<gmx::AwhHistory>();
             awhHistory      = awhHistoryLocal.get();
         }
 
@@ -1672,39 +1917,34 @@ static int do_cpt_awh(XDR *xd, gmx_bool bRead,
 }
 
 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)
             {
@@ -1714,14 +1954,13 @@ static int do_cpt_files(XDR *xd, gmx_bool bRead,
             {
                 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;
@@ -1743,19 +1982,19 @@ static int do_cpt_files(XDR *xd, gmx_bool bRead,
         }
         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.checksumSize,
                            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", outputfile.checksum.size(), outputfile.checksum.data(), list) != 0)
             {
                 return -1;
             }
         }
         else
         {
-            outputfiles[i].chksum_size = -1;
+            outputfile.checksumSize = -1;
         }
     }
     return 0;
@@ -1763,28 +2002,17 @@ static int do_cpt_files(XDR *xd, gmx_bool bRead,
 
 
 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,
-                      gmx_int64_t step, double t,
+                      int64_t step, double t,
                       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))
@@ -1811,16 +2039,16 @@ void write_checkpoint(const char *fn, gmx_bool bNumberAndKeep,
     snew(fntemp, std::strlen(fn));
     std::strcpy(fntemp, fn);
 #endif
-    gmx_format_current_time(timebuf, STRLEN);
+    std::string timebuf = gmx_format_current_time();
 
     if (fplog)
     {
         fprintf(fplog, "Writing checkpoint, step %s at %s\n\n",
-                gmx_step_str(step, buf), timebuf);
+                gmx_step_str(step, buf), timebuf.c_str());
     }
 
     /* 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");
 
@@ -1860,6 +2088,14 @@ void write_checkpoint(const char *fn, gmx_bool bNumberAndKeep,
         }
     }
 
+    PullHistory *pullHist         = observablesHistory->pullHistory.get();
+    int          flagsPullHistory = 0;
+    if (pullHist != nullptr && (pullHist->numValuesInXSum > 0 || pullHist->numValuesInFSum > 0))
+    {
+        flagsPullHistory |= (1<<epullhPULL_NUMCOORDINATES);
+        flagsPullHistory |= ((1<<epullhPULL_NUMGROUPS) | (1<<epullhPULL_NUMVALUESINXSUM) | (1<<epullhPULL_NUMVALUESINFSUM));
+    }
+
     int flags_dfh;
     if (bExpanded)
     {
@@ -1881,7 +2117,7 @@ void write_checkpoint(const char *fn, gmx_bool bNumberAndKeep,
     }
 
     int flags_awhh = 0;
-    if (state->awhHistory != nullptr && state->awhHistory->bias.size() > 0)
+    if (state->awhHistory != nullptr && !state->awhHistory->bias.empty())
     {
         flags_awhh |= ( (1 << eawhhIN_INITIAL) |
                         (1 << eawhhEQUILIBRATEHISTOGRAM) |
@@ -1900,53 +2136,49 @@ void write_checkpoint(const char *fn, gmx_bool bNumberAndKeep,
      * 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);
 
-    int             nlambda     = (state->dfhist ? state->dfhist->nlambda : 0);
+    edsamhistory_t          *edsamhist   = observablesHistory->edsamHistory.get();
+    int                      nED         = (edsamhist ? edsamhist->nED : 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);
 
-    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);
+    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,
+        flagsPullHistory, flags_dfh, flags_awhh,
+        nED, eSwapCoords
+    };
+    std::strcpy(headerContents.version, gmx_version());
+    std::strcpy(headerContents.fprog, gmx::getProgramContext().fullBinaryPath());
+    std::strcpy(headerContents.ftime, timebuf.c_str());
+    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_enerhist(gmx_fio_getxdr(fp), FALSE, flags_enh, enerhist, nullptr) < 0)  ||
+        (doCptPullHist(gmx_fio_getxdr(fp), FALSE, flagsPullHistory, pullHist, StatePart::pullHistory, nullptr) < 0)  ||
         (do_cpt_df_hist(gmx_fio_getxdr(fp), flags_dfh, nlambda, &state->dfhist, 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_awh(gmx_fio_getxdr(fp), FALSE, flags_awhh, state->awhHistory.get(), nullptr) < 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
@@ -1967,7 +2199,7 @@ void write_checkpoint(const char *fn, gmx_bool bNumberAndKeep,
         }
         else
         {
-            gmx_warning(buf);
+            gmx_warning("%s", buf);
         }
     }
 
@@ -1988,19 +2220,22 @@ void write_checkpoint(const char *fn, gmx_bool bNumberAndKeep,
             buf[std::strlen(fn) - std::strlen(ftp2ext(fn2ftp(fn))) - 1] = '\0';
             std::strcat(buf, "_prev");
             std::strcat(buf, fn+std::strlen(fn) - std::strlen(ftp2ext(fn2ftp(fn))) - 1);
-#ifndef GMX_FAHCORE
-            /* we copy here so that if something goes wrong between now and
-             * the rename below, there's always a state.cpt.
-             * If renames are atomic (such as in POSIX systems),
-             * this copying should be unneccesary.
-             */
-            gmx_file_copy(fn, buf, FALSE);
-            /* We don't really care if this fails:
-             * there's already a new checkpoint.
-             */
-#else
-            gmx_file_rename(fn, buf);
-#endif
+            if (!GMX_FAHCORE)
+            {
+                /* we copy here so that if something goes wrong between now and
+                 * the rename below, there's always a state.cpt.
+                 * If renames are atomic (such as in POSIX systems),
+                 * this copying should be unneccesary.
+                 */
+                gmx_file_copy(fn, buf, FALSE);
+                /* We don't really care if this fails:
+                 * there's already a new checkpoint.
+                 */
+            }
+            else
+            {
+                gmx_file_rename(fn, buf);
+            }
         }
         if (gmx_file_rename(fntemp, fn) != 0)
         {
@@ -2009,10 +2244,9 @@ void write_checkpoint(const char *fn, gmx_bool bNumberAndKeep,
     }
 #endif  /* GMX_NO_RENAME */
 
-    sfree(outputfiles);
     sfree(fntemp);
 
-#ifdef GMX_FAHCORE
+#if GMX_FAHCORE
     /*code for alternate checkpointing scheme.  moved from top of loop over
        steps */
     fcRequestCheckPoint();
@@ -2023,60 +2257,43 @@ void write_checkpoint(const char *fn, gmx_bool bNumberAndKeep,
 #endif /* end GMX_FAHCORE block */
 }
 
-static void print_flag_mismatch(FILE *fplog, int sflags, int fflags)
+static void check_int(FILE *fplog, const char *type, int p, int f, gmx_bool *mm)
 {
-    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++)
+    bool foundMismatch = (p != f);
+    if (!foundMismatch)
     {
-        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");
-        }
+        return;
     }
-}
-
-static void check_int(FILE *fplog, const char *type, int p, int f, gmx_bool *mm)
-{
-    FILE *fp = fplog ? fplog : stderr;
-
-    if (p != f)
+    *mm = TRUE;
+    if (fplog)
     {
-        fprintf(fp, "  %s mismatch,\n", type);
-        fprintf(fp, "    current program: %d\n", p);
-        fprintf(fp, "    checkpoint file: %d\n", f);
-        fprintf(fp, "\n");
-        *mm = TRUE;
+        fprintf(fplog, "  %s mismatch,\n", type);
+        fprintf(fplog, "    current program: %d\n", p);
+        fprintf(fplog, "    checkpoint file: %d\n", f);
+        fprintf(fplog, "\n");
     }
 }
 
 static void check_string(FILE *fplog, const char *type, const char *p,
                          const char *f, gmx_bool *mm)
 {
-    FILE *fp = fplog ? fplog : stderr;
-
-    if (std::strcmp(p, f) != 0)
+    bool foundMismatch = (std::strcmp(p, f) != 0);
+    if (!foundMismatch)
+    {
+        return;
+    }
+    *mm = TRUE;
+    if (fplog)
     {
-        fprintf(fp, "  %s mismatch,\n", type);
-        fprintf(fp, "    current program: %s\n", p);
-        fprintf(fp, "    checkpoint file: %s\n", f);
-        fprintf(fp, "\n");
-        *mm = TRUE;
+        fprintf(fplog, "  %s mismatch,\n", type);
+        fprintf(fplog, "    current program: %s\n", p);
+        fprintf(fplog, "    checkpoint file: %s\n", f);
+        fprintf(fplog, "\n");
     }
 }
 
-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
@@ -2084,16 +2301,15 @@ static void check_match(FILE *fplog,
      * 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[] =
             "You are continuing a simulation with a different precision. Not matching\n"
             "single/double precision will lead to precision or performance loss.\n";
-        fprintf(stderr, "%s\n", msg_precision_difference);
         if (fplog)
         {
             fprintf(fplog, "%s\n", msg_precision_difference);
@@ -2104,28 +2320,30 @@ static void check_match(FILE *fplog,
 
     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, "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);
         }
     }
 
@@ -2138,7 +2356,7 @@ static void check_match(FILE *fplog,
         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[] =
@@ -2154,13 +2372,8 @@ static void check_match(FILE *fplog,
             "GROMACS patchlevel, binary or parallel settings differ from previous run.\n"
             "Continuation is exact, but not guaranteed to be binary identical.\n";
 
-        const char msg_logdetails[] =
-            "See the log file for details.\n";
-
         if (majorVersionDiffers)
         {
-            fprintf(stderr, "%s%s\n", msg_major_version_difference, fplog ? msg_logdetails : "");
-
             if (fplog)
             {
                 fprintf(fplog, "%s\n", msg_major_version_difference);
@@ -2169,7 +2382,6 @@ static void check_match(FILE *fplog,
         else if (reproducibilityRequested)
         {
             /* Major & minor versions match at least, but something is different. */
-            fprintf(stderr, "%s%s\n", msg_mismatch_notice, fplog ? msg_logdetails : "");
             if (fplog)
             {
                 fprintf(fplog, "%s\n", msg_mismatch_notice);
@@ -2178,250 +2390,269 @@ static void check_match(FILE *fplog,
     }
 }
 
-static void read_checkpoint(const char *fn, FILE **pfplog,
+static void lockLogFile(FILE       *fplog,
+                        t_fileio   *chksum_file,
+                        const char *filename,
+                        bool        bForceAppend)
+{
+    /* Note that there are systems where the lock operation
+     * will succeed, but a second process can also lock the file.
+     * We should probably try to detect this.
+     */
+#if defined __native_client__
+    errno = ENOSYS;
+    if (true)
+#elif GMX_NATIVE_WINDOWS
+    if (_locking(fileno(gmx_fio_getfp(chksum_file)), _LK_NBLCK, LONG_MAX) == -1)
+#else
+    struct flock         fl; /* don't initialize here: the struct order is OS
+                                dependent! */
+    fl.l_type   = F_WRLCK;
+    fl.l_whence = SEEK_SET;
+    fl.l_start  = 0;
+    fl.l_len    = 0;
+    fl.l_pid    = 0;
+
+    if (fcntl(fileno(gmx_fio_getfp(chksum_file)), F_SETLK, &fl) == -1)
+#endif
+    {
+        if (errno == ENOSYS)
+        {
+            if (!bForceAppend)
+            {
+                gmx_fatal(FARGS, "File locking is not supported on this system. Use -noappend or specify -append explicitly to append anyhow.");
+            }
+            else
+            {
+                if (fplog)
+                {
+                    fprintf(fplog, "\nNOTE: File locking not supported on this system, will not lock %s\n\n", filename);
+                }
+            }
+        }
+        else if (errno == EACCES || errno == EAGAIN)
+        {
+            gmx_fatal(FARGS, "Failed to lock: %s. Already running "
+                      "simulation?", filename);
+        }
+        else
+        {
+            gmx_fatal(FARGS, "Failed to lock: %s. %s.",
+                      filename, std::strerror(errno));
+        }
+    }
+}
+
+//! Check whether chksum_file output file has a checksum that matches \c outputfile from the checkpoint.
+static void checkOutputFile(t_fileio                  *chksum_file,
+                            const gmx_file_position_t &outputfile)
+{
+    /* compute md5 chksum */
+    std::array<unsigned char, 16> digest;
+    if (outputfile.checksumSize != -1)
+    {
+        if (gmx_fio_get_file_md5(chksum_file, outputfile.offset,
+                                 &digest) != outputfile.checksumSize) /*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.",
+                      outputfile.checksumSize,
+                      outputfile.filename);
+        }
+    }
+
+    /* compare md5 chksum */
+    if (outputfile.checksumSize != -1 &&
+        digest != outputfile.checksum)
+    {
+        if (debug)
+        {
+            fprintf(debug, "chksum for %s: ", outputfile.filename);
+            for (int 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.",
+                  outputfile.filename);
+    }
+}
+
+static void read_checkpoint(const char *fn, t_fileio *logfio,
                             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                  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];
-#if !defined __native_client__ && !GMX_NATIVE_WINDOWS
-    struct flock         fl; /* don't initialize here: the struct order is OS
-                                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;
-    fl.l_start  = 0;
-    fl.l_len    = 0;
-    fl.l_pid    = 0;
-#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);
-    }
-
-    /* This will not be written if we do appending, since fplog is still NULL then */
+    // If we are appending, then we don't write to the open log file
+    // because we still need to compute a checksum for it. Otherwise
+    // we report to the new log file about the checkpoint file that we
+    // are reading from.
+    FILE *fplog = bAppendOutputFiles ? nullptr : gmx_fio_getfp(logfio);
     if (fplog)
     {
         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 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)) != 0) ||
+                  ((headerContents->flags_eks & (1<<eeksEKINF)) != 0) ||
+                  ((headerContents->flags_eks & (1<<eeksEKINO)) != 0) ||
+                  (((headerContents->flags_eks & (1<<eeksEKINSCALEF)) |
+                    (headerContents->flags_eks & (1<<eeksEKINSCALEH)) |
+                    (headerContents->flags_eks & (1<<eeksVSCALE))) != 0));
 
-    if (flags_enh && observablesHistory->energyHistory == nullptr)
+    if (headerContents->flags_enh && observablesHistory->energyHistory == nullptr)
     {
-        observablesHistory->energyHistory = std::unique_ptr<energyhistory_t>(new energyhistory_t {});
+        observablesHistory->energyHistory = gmx::compat::make_unique<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->flagsPullHistory)
     {
-        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(stderr, "\nWARNING: %s\n\n", warn);
-        if (fplog)
+        if (observablesHistory->pullHistory == nullptr)
         {
-            fprintf(fplog, "\nWARNING: %s\n\n", warn);
+            observablesHistory->pullHistory = gmx::compat::make_unique<PullHistory>();
         }
-        observablesHistory->energyHistory->nsum     = *step;
-        observablesHistory->energyHistory->nsum_sim = *step;
+        ret = doCptPullHist(gmx_fio_getxdr(fp), TRUE,
+                            headerContents->flagsPullHistory, observablesHistory->pullHistory.get(), StatePart::pullHistory, nullptr);
+        if (ret)
+        {
+            cp_error();
+        }
+    }
+
+    if (headerContents->file_version < 6)
+    {
+        gmx_fatal(FARGS, "Continuing from checkpoint files written before GROMACS 4.5 is not supported");
     }
 
-    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 {});
+        observablesHistory->edsamHistory = gmx::compat::make_unique<edsamhistory_t>(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());
+        state->awhHistory = std::make_shared<gmx::AwhHistory>();
     }
     ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE,
-                     flags_awhh, state->awhHistory.get(), NULL);
+                     headerContents->flags_awhh, state->awhHistory.get(), nullptr);
     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 {});
+        observablesHistory->swapHistory = gmx::compat::make_unique<swaphistory_t>(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();
@@ -2431,12 +2662,6 @@ static void read_checkpoint(const char *fn, FILE **pfplog,
         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
@@ -2446,6 +2671,10 @@ static void read_checkpoint(const char *fn, FILE **pfplog,
      */
     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
@@ -2454,134 +2683,59 @@ static void read_checkpoint(const char *fn, FILE **pfplog,
             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++)
+        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");
-
-#else
-            chksum_file = gmx_fio_open(outputfiles[i].filename, "r+");
+        }
 
-            /* lock log file */
-            if (i == 0)
+        // Handle the log file separately - it comes first in the
+        // list, so we make sure that we retain a lock on the open
+        // file that is never lifted after the checksum is calculated.
+        {
+            const gmx_file_position_t &logOutputFile   = outputfiles[0];
+            if (!GMX_FAHCORE)
             {
-                /* Note that there are systems where the lock operation
-                 * will succeed, but a second process can also lock the file.
-                 * We should probably try to detect this.
-                 */
-#if defined __native_client__
-                errno = ENOSYS;
-                if (1)
+                lockLogFile(fplog, logfio, logOutputFile.filename, bForceAppend);
+                checkOutputFile(logfio, logOutputFile);
 
-#elif GMX_NATIVE_WINDOWS
-                if (_locking(fileno(gmx_fio_getfp(chksum_file)), _LK_NBLCK, LONG_MAX) == -1)
-#else
-                if (fcntl(fileno(gmx_fio_getfp(chksum_file)), F_SETLK, &fl) == -1)
-#endif
-                {
-                    if (errno == ENOSYS)
-                    {
-                        if (!bForceAppend)
-                        {
-                            gmx_fatal(FARGS, "File locking is not supported on this system. Use -noappend or specify -append explicitly to append anyhow.");
-                        }
-                        else
-                        {
-                            fprintf(stderr, "\nNOTE: File locking is not supported on this system, will not lock %s\n\n", outputfiles[i].filename);
-                            if (fplog)
-                            {
-                                fprintf(fplog, "\nNOTE: File locking not supported on this system, will not lock %s\n\n", outputfiles[i].filename);
-                            }
-                        }
-                    }
-                    else if (errno == EACCES || errno == EAGAIN)
-                    {
-                        gmx_fatal(FARGS, "Failed to lock: %s. Already running "
-                                  "simulation?", outputfiles[i].filename);
-                    }
-                    else
-                    {
-                        gmx_fatal(FARGS, "Failed to lock: %s. %s.",
-                                  outputfiles[i].filename, std::strerror(errno));
-                    }
-                }
-            }
-
-            /* compute md5 chksum */
-            if (outputfiles[i].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*/
-                {
-                    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);
-                }
-            }
-            if (i == 0)  /*log file needs to be seeked in case we need to truncate (other files are truncated below)*/
-            {
-                if (gmx_fio_seek(chksum_file, outputfiles[i].offset))
+                if (gmx_fio_seek(logfio, logOutputFile.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 */
-            {
-                *pfplog = gmx_fio_getfp(chksum_file);
-            }
-            else
+        }
+        if (!GMX_FAHCORE)
+        {
+            // Now handle the remaining outputfiles
+            gmx::ArrayRef<gmx_file_position_t> outputfilesView(outputfiles);
+            for (const auto &outputfile : outputfilesView.subArray(1, outputfiles.size()-1))
             {
+                t_fileio *chksum_file = gmx_fio_open(outputfile.filename, "r+");
+                checkOutputFile(chksum_file, outputfile);
                 gmx_fio_close(chksum_file);
-            }
-#ifndef GMX_FAHCORE
-            /* compare md5 chksum */
-            if (outputfiles[i].chksum_size != -1 &&
-                memcmp(digest, outputfiles[i].chksum, 16) != 0)
-            {
-                if (debug)
+
+                if (!GMX_NATIVE_WINDOWS)
                 {
-                    fprintf(debug, "chksum for %s: ", outputfiles[i].filename);
-                    for (j = 0; j < 16; j++)
+                    /* For FAHCORE, we do this elsewhere*/
+                    rc = gmx_truncate(outputfile.filename, outputfile.offset);
+                    if (rc != 0)
                     {
-                        fprintf(debug, "%02x", digest[j]);
+                        gmx_fatal(FARGS, "Truncation of file %s failed. Cannot do appending because of this failure.", outputfile.filename);
                     }
-                    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);
-            }
-#endif
-
-
-            if (i != 0) /*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);
-                if (rc != 0)
-                {
-                    gmx_fatal(FARGS, "Truncation of file %s failed. Cannot do appending because of this failure.", outputfiles[i].filename);
-                }
-#endif
             }
         }
     }
-
-    sfree(outputfiles);
 }
 
 
-void load_checkpoint(const char *fn, FILE **fplog,
+void load_checkpoint(const char *fn, t_fileio *logfio,
                      const t_commrec *cr, const ivec dd_nc,
                      t_inputrec *ir, t_state *state,
                      gmx_bool *bReadEkin,
@@ -2589,22 +2743,21 @@ void load_checkpoint(const char *fn, FILE **fplog,
                      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,
+        read_checkpoint(fn, logfio,
                         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;
@@ -2612,13 +2765,13 @@ void load_checkpoint(const char *fn, FILE **fplog,
     // pass a checkpoint written by an normal completion to a restart,
     // mdrun will read all input, does some work but no steps, and
     // write successful output. But perhaps that is not desirable.
-    if ((ir->nsteps >= 0) && (ir->nsteps < step))
+    if ((ir->nsteps >= 0) && (ir->nsteps < headerContents.step))
     {
         // Note that we do not intend to support the use of mdrun
         // -nsteps to circumvent this condition.
         char nstepsString[STEPSTRSIZE], stepString[STEPSTRSIZE];
         gmx_step_str(ir->nsteps, nstepsString);
-        gmx_step_str(step, stepString);
+        gmx_step_str(headerContents.step, stepString);
         gmx_fatal(FARGS, "The input requested %s steps, however the checkpoint "
                   "file has already reached step %s. The simulation will not "
                   "proceed, because either your simulation is already complete, "
@@ -2627,82 +2780,57 @@ void load_checkpoint(const char *fn, FILE **fplog,
     }
     if (ir->nsteps >= 0)
     {
-        ir->nsteps          += ir->init_step - step;
+        ir->nsteps          += ir->init_step - headerContents.step;
     }
-    ir->init_step        = step;
-    ir->simulation_part += 1;
+    ir->init_step        = headerContents.step;
+    ir->simulation_part  = headerContents.simulation_part + 1;
 }
 
 void read_checkpoint_part_and_step(const char  *filename,
                                    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;
+                                   int64_t     *step)
+{
     t_fileio *fp;
 
     if (filename == nullptr ||
         !gmx_fexist(filename) ||
-        (!(fp = gmx_fio_open(filename, "r"))))
+        ((fp = gmx_fio_open(filename, "r")) == nullptr))
     {
         *simulation_part = 0;
         *step            = 0;
         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 =
+                                 int64_t *step, double *t, t_state *state,
+                                 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();
@@ -2710,118 +2838,94 @@ static void read_checkpoint_data(t_fileio *fp, int *simulation_part,
 
     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);
+    PullHistory pullHist = {};
+    ret = doCptPullHist(gmx_fio_getxdr(fp), TRUE,
+                        headerContents.flagsPullHistory, &pullHist, StatePart::pullHistory, nullptr);
+    if (ret)
+    {
+        cp_error();
+    }
+
+    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(), nullptr);
     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;
+    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;
-    fr->step    = gmx_int64_to_int(step,
-                                   "conversion of checkpoint to trajectory");
+    fr->step    = int64_to_int(step,
+                               "conversion of checkpoint to trajectory");
     fr->bTime      = TRUE;
     fr->time       = t;
     fr->bLambda    = TRUE;
     fr->lambda     = state.lambda[efptFEP];
     fr->fep_state  = state.fep_state;
     fr->bAtoms     = FALSE;
-    fr->bX         = (state.flags & (1<<estX));
+    fr->bX         = ((state.flags & (1<<estX)) != 0);
     if (fr->bX)
     {
         fr->x   = makeRvecArray(state.x, state.natoms);
     }
-    fr->bV      = (state.flags & (1<<estV));
+    fr->bV      = ((state.flags & (1<<estV)) != 0);
     if (fr->bV)
     {
         fr->v   = makeRvecArray(state.v, state.natoms);
     }
     fr->bF      = FALSE;
-    fr->bBox    = (state.flags & (1<<estBOX));
+    fr->bBox    = ((state.flags & (1<<estBOX)) != 0);
     if (fr->bBox)
     {
         copy_mat(state.box, fr->box);
@@ -2831,36 +2935,24 @@ void read_checkpoint_trxframe(t_fileio *fp, t_trxframe *fr)
 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();
@@ -2868,40 +2960,48 @@ void list_checkpoint(const char *fn, FILE *out)
 
     energyhistory_t enerhist;
     ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE,
-                          flags_enh, &enerhist, out);
+                          headerContents.flags_enh, &enerhist, out);
+
+    if (ret == 0)
+    {
+        PullHistory pullHist = {};
+        ret = doCptPullHist(gmx_fio_getxdr(fp), TRUE,
+                            headerContents.flagsPullHistory, &pullHist, StatePart::pullHistory, 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)
@@ -2916,17 +3016,15 @@ void list_checkpoint(const char *fn, FILE *out)
 
 /* 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;
+    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?");