src/gromacs/fileio/checkpoint.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2008,2009,2010,2011,2012 by the GROMACS development team.
   5  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   6  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37
  38 /* The source code in this file should be thread-safe.
  39    Please keep it that way. */
  40 #include "gmxpre.h"
  41
  42 #include "checkpoint.h"
  43
  44 #include <cerrno>
  45 #include <cstdlib>
  46 #include <cstring>
  47
  48 #include <array>
  49 #include <memory>
  50
  51 #include "buildinfo.h"
  52 #include "gromacs/fileio/filetypes.h"
  53 #include "gromacs/fileio/gmxfio.h"
  54 #include "gromacs/fileio/gmxfio_xdr.h"
  55 #include "gromacs/fileio/xdr_datatype.h"
  56 #include "gromacs/fileio/xdrf.h"
  57 #include "gromacs/gmxlib/network.h"
  58 #include "gromacs/math/vec.h"
  59 #include "gromacs/math/vecdump.h"
  60 #include "gromacs/math/vectypes.h"
  61 #include "gromacs/mdtypes/awh_correlation_history.h"
  62 #include "gromacs/mdtypes/awh_history.h"
  63 #include "gromacs/mdtypes/checkpointdata.h"
  64 #include "gromacs/mdtypes/commrec.h"
  65 #include "gromacs/mdtypes/df_history.h"
  66 #include "gromacs/mdtypes/edsamhistory.h"
  67 #include "gromacs/mdtypes/energyhistory.h"
  68 #include "gromacs/mdtypes/inputrec.h"
  69 #include "gromacs/mdtypes/md_enums.h"
  70 #include "gromacs/mdtypes/mdrunoptions.h"
  71 #include "gromacs/mdtypes/observableshistory.h"
  72 #include "gromacs/mdtypes/pullhistory.h"
  73 #include "gromacs/mdtypes/state.h"
  74 #include "gromacs/mdtypes/swaphistory.h"
  75 #include "gromacs/modularsimulator/modularsimulator.h"
  76 #include "gromacs/trajectory/trajectoryframe.h"
  77 #include "gromacs/utility/arrayref.h"
  78 #include "gromacs/utility/baseversion.h"
  79 #include "gromacs/utility/cstringutil.h"
  80 #include "gromacs/utility/fatalerror.h"
  81 #include "gromacs/utility/futil.h"
  82 #include "gromacs/utility/gmxassert.h"
  83 #include "gromacs/utility/int64_to_int.h"
  84 #include "gromacs/utility/keyvaluetree.h"
  85 #include "gromacs/utility/keyvaluetreebuilder.h"
  86 #include "gromacs/utility/keyvaluetreeserializer.h"
  87 #include "gromacs/utility/mdmodulenotification.h"
  88 #include "gromacs/utility/programcontext.h"
  89 #include "gromacs/utility/smalloc.h"
  90 #include "gromacs/utility/sysinfo.h"
  91 #include "gromacs/utility/textwriter.h"
  92 #include "gromacs/utility/txtdump.h"
  93
  94 #define CPT_MAGIC1 171817
  95 #define CPT_MAGIC2 171819
  96
  97 namespace gmx
  98 {
  99
 100 template<typename ValueType>
 101 void readKvtCheckpointValue(compat::not_null<ValueType*> value,
 102                             const std::string&           name,
 103                             const std::string&           identifier,
 104                             const KeyValueTreeObject&    kvt)
 105 {
 106     const std::string key = identifier + "-" + name;
 107     if (!kvt.keyExists(key))
 108     {
 109         std::string errorMessage = "Cannot read requested checkpoint value " + key + " .";
 110         GMX_THROW(InternalError(errorMessage));
 111     }
 112     *value = kvt[key].cast<ValueType>();
 113 }
 114
 115 template void readKvtCheckpointValue(compat::not_null<std::int64_t*> value,
 116                                      const std::string&              name,
 117                                      const std::string&              identifier,
 118                                      const KeyValueTreeObject&       kvt);
 119 template void readKvtCheckpointValue(compat::not_null<real*>   value,
 120                                      const std::string&        name,
 121                                      const std::string&        identifier,
 122                                      const KeyValueTreeObject& kvt);
 123
 124 template<typename ValueType>
 125 void writeKvtCheckpointValue(const ValueType&          value,
 126                              const std::string&        name,
 127                              const std::string&        identifier,
 128                              KeyValueTreeObjectBuilder kvtBuilder)
 129 {
 130     kvtBuilder.addValue<ValueType>(identifier + "-" + name, value);
 131 }
 132
 133 template void writeKvtCheckpointValue(const std::int64_t&       value,
 134                                       const std::string&        name,
 135                                       const std::string&        identifier,
 136                                       KeyValueTreeObjectBuilder kvtBuilder);
 137 template void writeKvtCheckpointValue(const real&               value,
 138                                       const std::string&        name,
 139                                       const std::string&        identifier,
 140                                       KeyValueTreeObjectBuilder kvtBuilder);
 141
 142
 143 } // namespace gmx
 144
 145 /*! \brief Enum of values that describe the contents of a cpt file
 146  * whose format matches a version number
 147  *
 148  * The enum helps the code be more self-documenting and ensure merges
 149  * do not silently resolve when two patches make the same bump. When
 150  * adding new functionality, add a new element just above cptv_Count
 151  * in this enumeration, and write code below that does the right thing
 152  * according to the value of file_version.
 153  */
 154 enum cptv
 155 {
 156     cptv_Unknown = 17,                  /**< Version before numbering scheme */
 157     cptv_RemoveBuildMachineInformation, /**< remove functionality that makes mdrun builds non-reproducible */
 158     cptv_ComPrevStepAsPullGroupReference, /**< Allow using COM of previous step as pull group PBC reference */
 159     cptv_PullAverage,      /**< Added possibility to output average pull force and position */
 160     cptv_MdModules,        /**< Added checkpointing for MdModules */
 161     cptv_ModularSimulator, /**< Added checkpointing for modular simulator */
 162     cptv_Count             /**< the total number of cptv versions */
 163 };
 164
 165 /*! \brief Version number of the file format written to checkpoint
 166  * files by this version of the code.
 167  *
 168  * cpt_version should normally only be changed, via adding a new field
 169  * to cptv enumeration, when the header or footer format changes.
 170  *
 171  * The state data format itself is backward and forward compatible.
 172  * But old code can not read a new entry that is present in the file
 173  * (but can read a new format when new entries are not present).
 174  *
 175  * The cpt_version increases whenever the file format in the main
 176  * development branch changes, due to an extension of the cptv enum above.
 177  * Backward compatibility for reading old run input files is maintained
 178  * by checking this version number against that of the file and then using
 179  * the correct code path. */
 180 static const int cpt_version = cptv_Count - 1;
 181
 182
 183 const char* est_names[estNR] = { "FE-lambda",
 184                                  "box",
 185                                  "box-rel",
 186                                  "box-v",
 187                                  "pres_prev",
 188                                  "nosehoover-xi",
 189                                  "thermostat-integral",
 190                                  "x",
 191                                  "v",
 192                                  "sdx-unsupported",
 193                                  "CGp",
 194                                  "LD-rng-unsupported",
 195                                  "LD-rng-i-unsupported",
 196                                  "disre_initf",
 197                                  "disre_rm3tav",
 198                                  "orire_initf",
 199                                  "orire_Dtav",
 200                                  "svir_prev",
 201                                  "nosehoover-vxi",
 202                                  "v_eta",
 203                                  "vol0",
 204                                  "nhpres_xi",
 205                                  "nhpres_vxi",
 206                                  "fvir_prev",
 207                                  "fep_state",
 208                                  "MC-rng-unsupported",
 209                                  "MC-rng-i-unsupported",
 210                                  "barostat-integral" };
 211
 212 enum
 213 {
 214     eeksEKIN_N,
 215     eeksEKINH,
 216     eeksDEKINDL,
 217     eeksMVCOS,
 218     eeksEKINF,
 219     eeksEKINO,
 220     eeksEKINSCALEF,
 221     eeksEKINSCALEH,
 222     eeksVSCALE,
 223     eeksEKINTOTAL,
 224     eeksNR
 225 };
 226
 227 static const char* eeks_names[eeksNR] = { "Ekin_n",         "Ekinh",          "dEkindlambda",
 228                                           "mv_cos",         "Ekinf",          "Ekinh_old",
 229                                           "EkinScaleF_NHC", "EkinScaleH_NHC", "Vscale_NHC",
 230                                           "Ekin_Total" };
 231
 232 enum
 233 {
 234     eenhENERGY_N,
 235     eenhENERGY_AVER,
 236     eenhENERGY_SUM,
 237     eenhENERGY_NSUM,
 238     eenhENERGY_SUM_SIM,
 239     eenhENERGY_NSUM_SIM,
 240     eenhENERGY_NSTEPS,
 241     eenhENERGY_NSTEPS_SIM,
 242     eenhENERGY_DELTA_H_NN,
 243     eenhENERGY_DELTA_H_LIST,
 244     eenhENERGY_DELTA_H_STARTTIME,
 245     eenhENERGY_DELTA_H_STARTLAMBDA,
 246     eenhNR
 247 };
 248
 249 enum
 250 {
 251     epullhPULL_NUMCOORDINATES,
 252     epullhPULL_NUMGROUPS,
 253     epullhPULL_NUMVALUESINXSUM,
 254     epullhPULL_NUMVALUESINFSUM,
 255     epullhNR
 256 };
 257
 258 enum
 259 {
 260     epullcoordh_VALUE_REF_SUM,
 261     epullcoordh_VALUE_SUM,
 262     epullcoordh_DR01_SUM,
 263     epullcoordh_DR23_SUM,
 264     epullcoordh_DR45_SUM,
 265     epullcoordh_FSCAL_SUM,
 266     epullcoordh_DYNAX_SUM,
 267     epullcoordh_NR
 268 };
 269
 270 enum
 271 {
 272     epullgrouph_X_SUM,
 273     epullgrouph_NR
 274 };
 275
 276 static const char* eenh_names[eenhNR] = { "energy_n",
 277                                           "energy_aver",
 278                                           "energy_sum",
 279                                           "energy_nsum",
 280                                           "energy_sum_sim",
 281                                           "energy_nsum_sim",
 282                                           "energy_nsteps",
 283                                           "energy_nsteps_sim",
 284                                           "energy_delta_h_nn",
 285                                           "energy_delta_h_list",
 286                                           "energy_delta_h_start_time",
 287                                           "energy_delta_h_start_lambda" };
 288
 289 static const char* ePullhNames[epullhNR] = { "pullhistory_numcoordinates",
 290                                              "pullhistory_numgroups",
 291                                              "pullhistory_numvaluesinxsum",
 292                                              "pullhistory_numvaluesinfsum" };
 293
 294 /* free energy history variables -- need to be preserved over checkpoint */
 295 enum
 296 {
 297     edfhBEQUIL,
 298     edfhNATLAMBDA,
 299     edfhWLHISTO,
 300     edfhWLDELTA,
 301     edfhSUMWEIGHTS,
 302     edfhSUMDG,
 303     edfhSUMMINVAR,
 304     edfhSUMVAR,
 305     edfhACCUMP,
 306     edfhACCUMM,
 307     edfhACCUMP2,
 308     edfhACCUMM2,
 309     edfhTIJ,
 310     edfhTIJEMP,
 311     edfhNR
 312 };
 313 /* free energy history variable names  */
 314 static const char* edfh_names[edfhNR] = { "bEquilibrated",
 315                                           "N_at_state",
 316                                           "Wang-Landau Histogram",
 317                                           "Wang-Landau Delta",
 318                                           "Weights",
 319                                           "Free Energies",
 320                                           "minvar",
 321                                           "variance",
 322                                           "accumulated_plus",
 323                                           "accumulated_minus",
 324                                           "accumulated_plus_2",
 325                                           "accumulated_minus_2",
 326                                           "Tij",
 327                                           "Tij_empirical" };
 328
 329 /* AWH biasing history variables */
 330 enum
 331 {
 332     eawhhIN_INITIAL,
 333     eawhhEQUILIBRATEHISTOGRAM,
 334     eawhhHISTSIZE,
 335     eawhhNPOINTS,
 336     eawhhCOORDPOINT,
 337     eawhhUMBRELLAGRIDPOINT,
 338     eawhhUPDATELIST,
 339     eawhhLOGSCALEDSAMPLEWEIGHT,
 340     eawhhNUMUPDATES,
 341     eawhhFORCECORRELATIONGRID,
 342     eawhhNR
 343 };
 344
 345 static const char* eawhh_names[eawhhNR] = { "awh_in_initial", "awh_equilibrateHistogram",
 346                                             "awh_histsize",   "awh_npoints",
 347                                             "awh_coordpoint", "awh_umbrellaGridpoint",
 348                                             "awh_updatelist", "awh_logScaledSampleWeight",
 349                                             "awh_numupdates", "awh_forceCorrelationGrid" };
 350
 351 enum
 352 {
 353     epullsPREVSTEPCOM,
 354     epullsNR
 355 };
 356
 357 static const char* epull_prev_step_com_names[epullsNR] = { "Pull groups prev step COM" };
 358
 359
 360 //! Higher level vector element type, only used for formatting checkpoint dumps
 361 enum class CptElementType
 362 {
 363     integer,  //!< integer
 364     real,     //!< float or double, not linked to precision of type real
 365     real3,    //!< float[3] or double[3], not linked to precision of type real
 366     matrix3x3 //!< float[3][3] or double[3][3], not linked to precision of type real
 367 };
 368
 369 //! \brief Parts of the checkpoint state, only used for reporting
 370 enum class StatePart
 371 {
 372     microState,         //!< The microstate of the simulated system
 373     kineticEnergy,      //!< Kinetic energy, needed for T/P-coupling state
 374     energyHistory,      //!< Energy observable statistics
 375     freeEnergyHistory,  //!< Free-energy state and observable statistics
 376     accWeightHistogram, //!< Accelerated weight histogram method state
 377     pullState,          //!< COM of previous step.
 378     pullHistory         //!< Pull history statistics (sums since last written output)
 379 };
 380
 381 //! \brief Return the name of a checkpoint entry based on part and part entry
 382 static const char* entryName(StatePart part, int ecpt)
 383 {
 384     switch (part)
 385     {
 386         case StatePart::microState: return est_names[ecpt];
 387         case StatePart::kineticEnergy: return eeks_names[ecpt];
 388         case StatePart::energyHistory: return eenh_names[ecpt];
 389         case StatePart::freeEnergyHistory: return edfh_names[ecpt];
 390         case StatePart::accWeightHistogram: return eawhh_names[ecpt];
 391         case StatePart::pullState: return epull_prev_step_com_names[ecpt];
 392         case StatePart::pullHistory: return ePullhNames[ecpt];
 393     }
 394
 395     return nullptr;
 396 }
 397
 398 static void cp_warning(FILE* fp)
 399 {
 400     fprintf(fp, "\nWARNING: Checkpoint file is corrupted or truncated\n\n");
 401 }
 402
 403 [[noreturn]] static void cp_error()
 404 {
 405     gmx_fatal(FARGS, "Checkpoint file corrupted/truncated, or maybe you are out of disk space?");
 406 }
 407
 408 static void do_cpt_string_err(XDR* xd, const char* desc, gmx::ArrayRef<char> s, FILE* list)
 409 {
 410     char* data = s.data();
 411     if (xdr_string(xd, &data, s.size()) == 0)
 412     {
 413         cp_error();
 414     }
 415     if (list)
 416     {
 417         fprintf(list, "%s = %s\n", desc, data);
 418     }
 419 }
 420
 421 static int do_cpt_int(XDR* xd, const char* desc, int* i, FILE* list)
 422 {
 423     if (xdr_int(xd, i) == 0)
 424     {
 425         return -1;
 426     }
 427     if (list)
 428     {
 429         fprintf(list, "%s = %d\n", desc, *i);
 430     }
 431     return 0;
 432 }
 433
 434 static int do_cpt_u_chars(XDR* xd, const char* desc, int n, unsigned char* i, FILE* list)
 435 {
 436     if (list)
 437     {
 438         fprintf(list, "%s = ", desc);
 439     }
 440     bool_t res = 1;
 441     for (int j = 0; j < n && res; j++)
 442     {
 443         res &= xdr_u_char(xd, &i[j]);
 444         if (list)
 445         {
 446             fprintf(list, "%02x", i[j]);
 447         }
 448     }
 449     if (list)
 450     {
 451         fprintf(list, "\n");
 452     }
 453     if (res == 0)
 454     {
 455         return -1;
 456     }
 457
 458     return 0;
 459 }
 460
 461 static void do_cpt_int_err(XDR* xd, const char* desc, int* i, FILE* list)
 462 {
 463     if (do_cpt_int(xd, desc, i, list) < 0)
 464     {
 465         cp_error();
 466     }
 467 }
 468
 469 static void do_cpt_bool_err(XDR* xd, const char* desc, bool* b, FILE* list)
 470 {
 471     int i = static_cast<int>(*b);
 472
 473     if (do_cpt_int(xd, desc, &i, list) < 0)
 474     {
 475         cp_error();
 476     }
 477
 478     *b = (i != 0);
 479 }
 480
 481 static void do_cpt_step_err(XDR* xd, const char* desc, int64_t* i, FILE* list)
 482 {
 483     char buf[STEPSTRSIZE];
 484
 485     if (xdr_int64(xd, i) == 0)
 486     {
 487         cp_error();
 488     }
 489     if (list)
 490     {
 491         fprintf(list, "%s = %s\n", desc, gmx_step_str(*i, buf));
 492     }
 493 }
 494
 495 static void do_cpt_double_err(XDR* xd, const char* desc, double* f, FILE* list)
 496 {
 497     if (xdr_double(xd, f) == 0)
 498     {
 499         cp_error();
 500     }
 501     if (list)
 502     {
 503         fprintf(list, "%s = %f\n", desc, *f);
 504     }
 505 }
 506
 507 static void do_cpt_real_err(XDR* xd, real* f)
 508 {
 509 #if GMX_DOUBLE
 510     bool_t res = xdr_double(xd, f);
 511 #else
 512     bool_t res = xdr_float(xd, f);
 513 #endif
 514     if (res == 0)
 515     {
 516         cp_error();
 517     }
 518 }
 519
 520 static void do_cpt_n_rvecs_err(XDR* xd, const char* desc, int n, rvec f[], FILE* list)
 521 {
 522     for (int i = 0; i < n; i++)
 523     {
 524         for (int j = 0; j < DIM; j++)
 525         {
 526             do_cpt_real_err(xd, &f[i][j]);
 527         }
 528     }
 529
 530     if (list)
 531     {
 532         pr_rvecs(list, 0, desc, f, n);
 533     }
 534 }
 535
 536 template<typename T>
 537 struct xdr_type
 538 {
 539 };
 540
 541 template<>
 542 struct xdr_type<int>
 543 {
 544     static const int value = xdr_datatype_int;
 545 };
 546
 547 template<>
 548 struct xdr_type<float>
 549 {
 550     static const int value = xdr_datatype_float;
 551 };
 552
 553 template<>
 554 struct xdr_type<double>
 555 {
 556     static const int value = xdr_datatype_double;
 557 };
 558
 559 //! \brief Returns size in byte of an xdr_datatype
 560 static inline unsigned int sizeOfXdrType(int xdrType)
 561 {
 562     switch (xdrType)
 563     {
 564         case xdr_datatype_int: return sizeof(int);
 565         case xdr_datatype_float: return sizeof(float);
 566         case xdr_datatype_double: return sizeof(double);
 567         default: GMX_RELEASE_ASSERT(false, "XDR data type not implemented");
 568     }
 569
 570     return 0;
 571 }
 572
 573 //! \brief Returns the XDR process function for i/o of an XDR type
 574 static inline xdrproc_t xdrProc(int xdrType)
 575 {
 576     switch (xdrType)
 577     {
 578         case xdr_datatype_int: return reinterpret_cast<xdrproc_t>(xdr_int);
 579         case xdr_datatype_float: return reinterpret_cast<xdrproc_t>(xdr_float);
 580         case xdr_datatype_double: return reinterpret_cast<xdrproc_t>(xdr_double);
 581         default: GMX_RELEASE_ASSERT(false, "XDR data type not implemented");
 582     }
 583
 584     return nullptr;
 585 }
 586
 587 /*! \brief Lists or only reads an xdr vector from checkpoint file
 588  *
 589  * When list!=NULL reads and lists the \p nf vector elements of type \p xdrType.
 590  * The header for the print is set by \p part and \p ecpt.
 591  * The formatting of the printing is set by \p cptElementType.
 592  * When list==NULL only reads the elements.
 593  */
 594 static bool_t listXdrVector(XDR* xd, StatePart part, int ecpt, int nf, int xdrType, FILE* list, CptElementType cptElementType)
 595 {
 596     bool_t res = 0;
 597
 598     const unsigned int elemSize = sizeOfXdrType(xdrType);
 599     std::vector<char>  data(nf * elemSize);
 600     res = xdr_vector(xd, data.data(), nf, elemSize, xdrProc(xdrType));
 601
 602     if (list != nullptr)
 603     {
 604         switch (xdrType)
 605         {
 606             case xdr_datatype_int:
 607                 pr_ivec(list, 0, entryName(part, ecpt), reinterpret_cast<const int*>(data.data()), nf, TRUE);
 608                 break;
 609             case xdr_datatype_float:
 610 #if !GMX_DOUBLE
 611                 if (cptElementType == CptElementType::real3)
 612                 {
 613                     pr_rvecs(list, 0, entryName(part, ecpt), reinterpret_cast<const rvec*>(data.data()), nf / 3);
 614                 }
 615                 else
 616 #endif
 617                 {
 618                     /* Note: With double precision code dumping a single precision rvec will produce float iso rvec print, but that's a minor annoyance */
 619                     pr_fvec(list, 0, entryName(part, ecpt), reinterpret_cast<const float*>(data.data()), nf, TRUE);
 620                 }
 621                 break;
 622             case xdr_datatype_double:
 623 #if GMX_DOUBLE
 624                 if (cptElementType == CptElementType::real3)
 625                 {
 626                     pr_rvecs(list, 0, entryName(part, ecpt), reinterpret_cast<const rvec*>(data.data()), nf / 3);
 627                 }
 628                 else
 629 #endif
 630                 {
 631                     /* Note: With single precision code dumping a double precision rvec will produce float iso rvec print, but that's a minor annoyance */
 632                     pr_dvec(list, 0, entryName(part, ecpt), reinterpret_cast<const double*>(data.data()), nf, TRUE);
 633                 }
 634                 break;
 635             default: GMX_RELEASE_ASSERT(false, "Data type not implemented for listing");
 636         }
 637     }
 638
 639     return res;
 640 }
 641
 642 //! \brief Convert a double array, typed char*, to float
 643 gmx_unused static void convertArrayRealPrecision(const char* c, float* v, int n)
 644 {
 645     const double* d = reinterpret_cast<const double*>(c);
 646     for (int i = 0; i < n; i++)
 647     {
 648         v[i] = static_cast<float>(d[i]);
 649     }
 650 }
 651
 652 //! \brief Convert a float array, typed char*, to double
 653 static void convertArrayRealPrecision(const char* c, double* v, int n)
 654 {
 655     const float* f = reinterpret_cast<const float*>(c);
 656     for (int i = 0; i < n; i++)
 657     {
 658         v[i] = static_cast<double>(f[i]);
 659     }
 660 }
 661
 662 //! \brief Generate an error for trying to convert to integer
 663 static void convertArrayRealPrecision(const char gmx_unused* c, int gmx_unused* v, int gmx_unused n)
 664 {
 665     GMX_RELEASE_ASSERT(false,
 666                        "We only expect type mismatches between float and double, not integer");
 667 }
 668
 669 /*! \brief Low-level routine for reading/writing a vector of reals from/to file.
 670  *
 671  * This is the only routine that does the actually i/o of real vector,
 672  * all other routines are intermediate level routines for specific real
 673  * data types, calling this routine.
 674  * Currently this routine is (too) complex, since it handles both real *
 675  * and std::vector<real>. Using real * is deprecated and this routine
 676  * will simplify a lot when only std::vector needs to be supported.
 677  *
 678  * The routine is generic to vectors with different allocators,
 679  * because as part of reading a checkpoint there are vectors whose
 680  * size is not known until reading has progressed far enough, so a
 681  * resize method must be called.
 682  *
 683  * When not listing, we use either v or vector, depending on which is !=NULL.
 684  * If nval >= 0, nval is used; on read this should match the passed value.
 685  * If nval n<0, *nptr (with v) or vector->size() is used. On read using v,
 686  * the value is stored in nptr
 687  */
 688 template<typename T, typename AllocatorType>
 689 static int doVectorLow(XDR*                           xd,
 690                        StatePart                      part,
 691                        int                            ecpt,
 692                        int                            sflags,
 693                        int                            nval,
 694                        int*                           nptr,
 695                        T**                            v,
 696                        std::vector<T, AllocatorType>* vector,
 697                        FILE*                          list,
 698                        CptElementType                 cptElementType)
 699 {
 700     GMX_RELEASE_ASSERT(list != nullptr || (v != nullptr && vector == nullptr)
 701                                || (v == nullptr && vector != nullptr),
 702                        "Without list, we should have exactly one of v and vector != NULL");
 703
 704     bool_t res = 0;
 705
 706     int numElemInTheFile;
 707     if (list == nullptr)
 708     {
 709         if (nval >= 0)
 710         {
 711             GMX_RELEASE_ASSERT(nptr == nullptr, "With nval>=0 we should have nptr==NULL");
 712             numElemInTheFile = nval;
 713         }
 714         else
 715         {
 716             if (v != nullptr)
 717             {
 718                 GMX_RELEASE_ASSERT(nptr != nullptr, "With nval<0 we should have nptr!=NULL");
 719                 numElemInTheFile = *nptr;
 720             }
 721             else
 722             {
 723                 numElemInTheFile = vector->size();
 724             }
 725         }
 726     }
 727     /* Read/write the vector element count */
 728     res = xdr_int(xd, &numElemInTheFile);
 729     if (res == 0)
 730     {
 731         return -1;
 732     }
 733     /* Read/write the element data type */
 734     constexpr int xdrTypeInTheCode = xdr_type<T>::value;
 735     int           xdrTypeInTheFile = xdrTypeInTheCode;
 736     res                            = xdr_int(xd, &xdrTypeInTheFile);
 737     if (res == 0)
 738     {
 739         return -1;
 740     }
 741
 742     if (list == nullptr && (sflags & (1 << ecpt)))
 743     {
 744         if (nval >= 0)
 745         {
 746             if (numElemInTheFile != nval)
 747             {
 748                 gmx_fatal(FARGS,
 749                           "Count mismatch for state entry %s, code count is %d, file count is %d\n",
 750                           entryName(part, ecpt),
 751                           nval,
 752                           numElemInTheFile);
 753             }
 754         }
 755         else if (nptr != nullptr)
 756         {
 757             *nptr = numElemInTheFile;
 758         }
 759
 760         bool typesMatch = (xdrTypeInTheFile == xdrTypeInTheCode);
 761         if (!typesMatch)
 762         {
 763             char buf[STRLEN];
 764             sprintf(buf,
 765                     "mismatch for state entry %s, code precision is %s, file precision is %s",
 766                     entryName(part, ecpt),
 767                     xdr_datatype_names[xdrTypeInTheCode],
 768                     xdr_datatype_names[xdrTypeInTheFile]);
 769
 770             /* Matching int and real should never occur, but check anyhow */
 771             if (xdrTypeInTheFile == xdr_datatype_int || xdrTypeInTheCode == xdr_datatype_int)
 772             {
 773                 gmx_fatal(FARGS,
 774                           "Type %s: incompatible checkpoint formats or corrupted checkpoint file.",
 775                           buf);
 776             }
 777         }
 778
 779         T* vp;
 780         if (v != nullptr)
 781         {
 782             if (*v == nullptr)
 783             {
 784                 snew(*v, numElemInTheFile);
 785             }
 786             vp = *v;
 787         }
 788         else
 789         {
 790             /* This conditional ensures that we don't resize on write.
 791              * In particular in the state where this code was written
 792              * vector has a size of numElemInThefile and we
 793              * don't want to lose that padding here.
 794              */
 795             if (vector->size() < static_cast<unsigned int>(numElemInTheFile))
 796             {
 797                 vector->resize(numElemInTheFile);
 798             }
 799             vp = vector->data();
 800         }
 801
 802         char* vChar;
 803         if (typesMatch)
 804         {
 805             vChar = reinterpret_cast<char*>(vp);
 806         }
 807         else
 808         {
 809             snew(vChar, numElemInTheFile * sizeOfXdrType(xdrTypeInTheFile));
 810         }
 811         res = xdr_vector(
 812                 xd, vChar, numElemInTheFile, sizeOfXdrType(xdrTypeInTheFile), xdrProc(xdrTypeInTheFile));
 813         if (res == 0)
 814         {
 815             return -1;
 816         }
 817
 818         if (!typesMatch)
 819         {
 820             /* In the old code float-double conversion came for free.
 821              * In the new code we still support it, mainly because
 822              * the tip4p_continue regression test makes use of this.
 823              * It's an open question if we do or don't want to allow this.
 824              */
 825             convertArrayRealPrecision(vChar, vp, numElemInTheFile);
 826             sfree(vChar);
 827         }
 828     }
 829     else
 830     {
 831         res = listXdrVector(xd, part, ecpt, numElemInTheFile, xdrTypeInTheFile, list, cptElementType);
 832     }
 833
 834     return 0;
 835 }
 836
 837 //! \brief Read/Write a std::vector, on read checks the number of elements matches \p numElements, if specified.
 838 template<typename T>
 839 static int
 840 doVector(XDR* xd, StatePart part, int ecpt, int sflags, std::vector<T>* vector, FILE* list, int numElements = -1)
 841 {
 842     return doVectorLow<T>(
 843             xd, part, ecpt, sflags, numElements, nullptr, nullptr, vector, list, CptElementType::real);
 844 }
 845
 846 //! \brief Read/Write an ArrayRef<real>.
 847 static int doRealArrayRef(XDR* xd, StatePart part, int ecpt, int sflags, gmx::ArrayRef<real> vector, FILE* list)
 848 {
 849     real* v_real = vector.data();
 850     return doVectorLow<real, std::allocator<real>>(
 851             xd, part, ecpt, sflags, vector.size(), nullptr, &v_real, nullptr, list, CptElementType::real);
 852 }
 853
 854 //! Convert from view of RVec to view of real.
 855 static gmx::ArrayRef<real> realArrayRefFromRVecArrayRef(gmx::ArrayRef<gmx::RVec> ofRVecs)
 856 {
 857     return gmx::arrayRefFromArray<real>(reinterpret_cast<real*>(ofRVecs.data()), ofRVecs.size() * DIM);
 858 }
 859
 860 //! \brief Read/Write a PaddedVector whose value_type is RVec.
 861 template<typename PaddedVectorOfRVecType>
 862 static int
 863 doRvecVector(XDR* xd, StatePart part, int ecpt, int sflags, PaddedVectorOfRVecType* v, int numAtoms, FILE* list)
 864 {
 865     const int numReals = numAtoms * DIM;
 866
 867     if (list == nullptr)
 868     {
 869         GMX_RELEASE_ASSERT(
 870                 sflags & (1 << ecpt),
 871                 "When not listing, the flag for the entry should be set when requesting i/o");
 872         GMX_RELEASE_ASSERT(v->size() == numAtoms, "v should have sufficient size for numAtoms");
 873
 874         return doRealArrayRef(xd, part, ecpt, sflags, realArrayRefFromRVecArrayRef(makeArrayRef(*v)), list);
 875     }
 876     else
 877     {
 878         // Use the rebind facility to change the value_type of the
 879         // allocator from RVec to real.
 880         using realAllocator =
 881                 typename std::allocator_traits<typename PaddedVectorOfRVecType::allocator_type>::template rebind_alloc<real>;
 882         return doVectorLow<real, realAllocator>(
 883                 xd, part, ecpt, sflags, numReals, nullptr, nullptr, nullptr, list, CptElementType::real);
 884     }
 885 }
 886
 887 /* This function stores n along with the reals for reading,
 888  * but on reading it assumes that n matches the value in the checkpoint file,
 889  * a fatal error is generated when this is not the case.
 890  */
 891 static int do_cpte_reals(XDR* xd, StatePart part, int ecpt, int sflags, int n, real** v, FILE* list)
 892 {
 893     return doVectorLow<real, std::allocator<real>>(
 894             xd, part, ecpt, sflags, n, nullptr, v, nullptr, list, CptElementType::real);
 895 }
 896
 897 /* This function does the same as do_cpte_reals,
 898  * except that on reading it ignores the passed value of *n
 899  * and stores the value read from the checkpoint file in *n.
 900  */
 901 static int do_cpte_n_reals(XDR* xd, StatePart part, int ecpt, int sflags, int* n, real** v, FILE* list)
 902 {
 903     return doVectorLow<real, std::allocator<real>>(
 904             xd, part, ecpt, sflags, -1, n, v, nullptr, list, CptElementType::real);
 905 }
 906
 907 static int do_cpte_real(XDR* xd, StatePart part, int ecpt, int sflags, real* r, FILE* list)
 908 {
 909     return doVectorLow<real, std::allocator<real>>(
 910             xd, part, ecpt, sflags, 1, nullptr, &r, nullptr, list, CptElementType::real);
 911 }
 912
 913 static int do_cpte_ints(XDR* xd, StatePart part, int ecpt, int sflags, int n, int** v, FILE* list)
 914 {
 915     return doVectorLow<int, std::allocator<int>>(
 916             xd, part, ecpt, sflags, n, nullptr, v, nullptr, list, CptElementType::integer);
 917 }
 918
 919 static int do_cpte_int(XDR* xd, StatePart part, int ecpt, int sflags, int* i, FILE* list)
 920 {
 921     return do_cpte_ints(xd, part, ecpt, sflags, 1, &i, list);
 922 }
 923
 924 static int do_cpte_bool(XDR* xd, StatePart part, int ecpt, int sflags, bool* b, FILE* list)
 925 {
 926     int i   = static_cast<int>(*b);
 927     int ret = do_cpte_int(xd, part, ecpt, sflags, &i, list);
 928     *b      = (i != 0);
 929     return ret;
 930 }
 931
 932 static int do_cpte_doubles(XDR* xd, StatePart part, int ecpt, int sflags, int n, double** v, FILE* list)
 933 {
 934     return doVectorLow<double, std::allocator<double>>(
 935             xd, part, ecpt, sflags, n, nullptr, v, nullptr, list, CptElementType::real);
 936 }
 937
 938 static int do_cpte_double(XDR* xd, StatePart part, int ecpt, int sflags, double* r, FILE* list)
 939 {
 940     return do_cpte_doubles(xd, part, ecpt, sflags, 1, &r, list);
 941 }
 942
 943 static int do_cpte_matrix(XDR* xd, StatePart part, int ecpt, int sflags, matrix v, FILE* list)
 944 {
 945     real* vr;
 946     int   ret;
 947
 948     vr  = &(v[0][0]);
 949     ret = doVectorLow<real, std::allocator<real>>(
 950             xd, part, ecpt, sflags, DIM * DIM, nullptr, &vr, nullptr, nullptr, CptElementType::matrix3x3);
 951
 952     if (list && ret == 0)
 953     {
 954         pr_rvecs(list, 0, entryName(part, ecpt), v, DIM);
 955     }
 956
 957     return ret;
 958 }
 959
 960
 961 static int do_cpte_nmatrix(XDR* xd, StatePart part, int ecpt, int sflags, int n, real** v, FILE* list)
 962 {
 963     int  i;
 964     int  ret, reti;
 965     char name[CPTSTRLEN];
 966
 967     ret = 0;
 968     if (v == nullptr)
 969     {
 970         snew(v, n);
 971     }
 972     for (i = 0; i < n; i++)
 973     {
 974         reti = doVectorLow<real, std::allocator<real>>(
 975                 xd, part, ecpt, sflags, n, nullptr, &(v[i]), nullptr, nullptr, CptElementType::matrix3x3);
 976         if (list && reti == 0)
 977         {
 978             sprintf(name, "%s[%d]", entryName(part, ecpt), i);
 979             pr_reals(list, 0, name, v[i], n);
 980         }
 981         if (reti != 0)
 982         {
 983             ret = reti;
 984         }
 985     }
 986     return ret;
 987 }
 988
 989 static int do_cpte_matrices(XDR* xd, StatePart part, int ecpt, int sflags, int n, matrix** v, FILE* list)
 990 {
 991     bool_t  res = 0;
 992     matrix *vp, *va = nullptr;
 993     real*   vr;
 994     int     nf, i, j, k;
 995     int     ret;
 996
 997     nf  = n;
 998     res = xdr_int(xd, &nf);
 999     if (res == 0)
1000     {
1001         return -1;
1002     }
1003     if (list == nullptr && nf != n)
1004     {
1005         gmx_fatal(FARGS,
1006                   "Count mismatch for state entry %s, code count is %d, file count is %d\n",
1007                   entryName(part, ecpt),
1008                   n,
1009                   nf);
1010     }
1011     if (list || !(sflags & (1 << ecpt)))
1012     {
1013         snew(va, nf);
1014         vp = va;
1015     }
1016     else
1017     {
1018         if (*v == nullptr)
1019         {
1020             snew(*v, nf);
1021         }
1022         vp = *v;
1023     }
1024     snew(vr, nf * DIM * DIM);
1025     for (i = 0; i < nf; i++)
1026     {
1027         for (j = 0; j < DIM; j++)
1028         {
1029             for (k = 0; k < DIM; k++)
1030             {
1031                 vr[(i * DIM + j) * DIM + k] = vp[i][j][k];
1032             }
1033         }
1034     }
1035     ret = doVectorLow<real, std::allocator<real>>(
1036             xd, part, ecpt, sflags, nf * DIM * DIM, nullptr, &vr, nullptr, nullptr, CptElementType::matrix3x3);
1037     for (i = 0; i < nf; i++)
1038     {
1039         for (j = 0; j < DIM; j++)
1040         {
1041             for (k = 0; k < DIM; k++)
1042             {
1043                 vp[i][j][k] = vr[(i * DIM + j) * DIM + k];
1044             }
1045         }
1046     }
1047     sfree(vr);
1048
1049     if (list && ret == 0)
1050     {
1051         for (i = 0; i < nf; i++)
1052         {
1053             pr_rvecs(list, 0, entryName(part, ecpt), vp[i], DIM);
1054         }
1055     }
1056     if (va)
1057     {
1058         sfree(va);
1059     }
1060
1061     return ret;
1062 }
1063
1064 static void do_cpt_header(XDR* xd, gmx_bool bRead, FILE* list, CheckpointHeaderContents* contents)
1065 {
1066     bool_t res = 0;
1067     int    magic;
1068
1069     if (bRead)
1070     {
1071         magic = -1;
1072     }
1073     else
1074     {
1075         magic = CPT_MAGIC1;
1076     }
1077     res = xdr_int(xd, &magic);
1078     if (res == 0)
1079     {
1080         gmx_fatal(FARGS,
1081                   "The checkpoint file is empty/corrupted, or maybe you are out of disk space?");
1082     }
1083     if (magic != CPT_MAGIC1)
1084     {
1085         gmx_fatal(FARGS,
1086                   "Start of file magic number mismatch, checkpoint file has %d, should be %d\n"
1087                   "The checkpoint file is corrupted or not a checkpoint file",
1088                   magic,
1089                   CPT_MAGIC1);
1090     }
1091     char fhost[255];
1092     if (!bRead)
1093     {
1094         gmx_gethostname(fhost, 255);
1095     }
1096     do_cpt_string_err(xd, "GROMACS version", contents->version, list);
1097     // The following fields are no longer ever written with meaningful
1098     // content, but because they precede the file version, there is no
1099     // good way for new code to read the old and new formats, nor a
1100     // good way for old code to avoid giving an error while reading a
1101     // new format. So we read and write a field that no longer has a
1102     // purpose.
1103     do_cpt_string_err(xd, "GROMACS build time UNUSED", contents->btime_UNUSED, list);
1104     do_cpt_string_err(xd, "GROMACS build user UNUSED", contents->buser_UNUSED, list);
1105     do_cpt_string_err(xd, "GROMACS build host UNUSED", contents->bhost_UNUSED, list);
1106     do_cpt_string_err(xd, "generating program", contents->fprog, list);
1107     do_cpt_string_err(xd, "generation time", contents->ftime, list);
1108     contents->file_version = cpt_version;
1109     do_cpt_int_err(xd, "checkpoint file version", &contents->file_version, list);
1110     if (contents->file_version > cpt_version)
1111     {
1112         gmx_fatal(FARGS,
1113                   "Attempting to read a checkpoint file of version %d with code of version %d\n",
1114                   contents->file_version,
1115                   cpt_version);
1116     }
1117     if (contents->file_version >= 13)
1118     {
1119         do_cpt_int_err(xd, "GROMACS double precision", &contents->double_prec, list);
1120     }
1121     else
1122     {
1123         contents->double_prec = -1;
1124     }
1125     if (contents->file_version >= 12)
1126     {
1127         do_cpt_string_err(xd, "generating host", fhost, list);
1128     }
1129     do_cpt_int_err(xd, "#atoms", &contents->natoms, list);
1130     do_cpt_int_err(xd, "#T-coupling groups", &contents->ngtc, list);
1131     if (contents->file_version >= 10)
1132     {
1133         do_cpt_int_err(xd, "#Nose-Hoover T-chains", &contents->nhchainlength, list);
1134     }
1135     else
1136     {
1137         contents->nhchainlength = 1;
1138     }
1139     if (contents->file_version >= 11)
1140     {
1141         do_cpt_int_err(xd, "#Nose-Hoover T-chains for barostat ", &contents->nnhpres, list);
1142     }
1143     else
1144     {
1145         contents->nnhpres = 0;
1146     }
1147     if (contents->file_version >= 14)
1148     {
1149         do_cpt_int_err(xd, "# of total lambda states ", &contents->nlambda, list);
1150     }
1151     else
1152     {
1153         contents->nlambda = 0;
1154     }
1155     do_cpt_int_err(xd, "integrator", &contents->eIntegrator, list);
1156     if (contents->file_version >= 3)
1157     {
1158         do_cpt_int_err(xd, "simulation part #", &contents->simulation_part, list);
1159     }
1160     else
1161     {
1162         contents->simulation_part = 1;
1163     }
1164     if (contents->file_version >= 5)
1165     {
1166         do_cpt_step_err(xd, "step", &contents->step, list);
1167     }
1168     else
1169     {
1170         int idum = 0;
1171         do_cpt_int_err(xd, "step", &idum, list);
1172         contents->step = static_cast<int64_t>(idum);
1173     }
1174     do_cpt_double_err(xd, "t", &contents->t, list);
1175     do_cpt_int_err(xd, "#PP-ranks", &contents->nnodes, list);
1176     do_cpt_int_err(xd, "dd_nc[x]", &contents->dd_nc[XX], list);
1177     do_cpt_int_err(xd, "dd_nc[y]", &contents->dd_nc[YY], list);
1178     do_cpt_int_err(xd, "dd_nc[z]", &contents->dd_nc[ZZ], list);
1179     do_cpt_int_err(xd, "#PME-only ranks", &contents->npme, list);
1180     do_cpt_int_err(xd, "state flags", &contents->flags_state, list);
1181     if (contents->file_version >= 4)
1182     {
1183         do_cpt_int_err(xd, "ekin data flags", &contents->flags_eks, list);
1184         do_cpt_int_err(xd, "energy history flags", &contents->flags_enh, list);
1185     }
1186     else
1187     {
1188         contents->flags_eks   = 0;
1189         contents->flags_enh   = (contents->flags_state >> (estORIRE_DTAV + 1));
1190         contents->flags_state = (contents->flags_state
1191                                  & ~((1 << (estORIRE_DTAV + 1)) | (1 << (estORIRE_DTAV + 2))
1192                                      | (1 << (estORIRE_DTAV + 3))));
1193     }
1194     if (contents->file_version >= 14)
1195     {
1196         do_cpt_int_err(xd, "df history flags", &contents->flags_dfh, list);
1197     }
1198     else
1199     {
1200         contents->flags_dfh = 0;
1201     }
1202
1203     if (contents->file_version >= 15)
1204     {
1205         do_cpt_int_err(xd, "ED data sets", &contents->nED, list);
1206     }
1207     else
1208     {
1209         contents->nED = 0;
1210     }
1211
1212     if (contents->file_version >= 16)
1213     {
1214         do_cpt_int_err(xd, "swap", &contents->eSwapCoords, list);
1215     }
1216     else
1217     {
1218         contents->eSwapCoords = eswapNO;
1219     }
1220
1221     if (contents->file_version >= 17)
1222     {
1223         do_cpt_int_err(xd, "AWH history flags", &contents->flags_awhh, list);
1224     }
1225     else
1226     {
1227         contents->flags_awhh = 0;
1228     }
1229
1230     if (contents->file_version >= 18)
1231     {
1232         do_cpt_int_err(xd, "pull history flags", &contents->flagsPullHistory, list);
1233     }
1234     else
1235     {
1236         contents->flagsPullHistory = 0;
1237     }
1238
1239     if (contents->file_version >= cptv_ModularSimulator)
1240     {
1241         do_cpt_bool_err(
1242                 xd, "Is modular simulator checkpoint", &contents->isModularSimulatorCheckpoint, list);
1243     }
1244     else
1245     {
1246         contents->isModularSimulatorCheckpoint = false;
1247     }
1248 }
1249
1250 static int do_cpt_footer(XDR* xd, int file_version)
1251 {
1252     bool_t res = 0;
1253     int    magic;
1254
1255     if (file_version >= 2)
1256     {
1257         magic = CPT_MAGIC2;
1258         res   = xdr_int(xd, &magic);
1259         if (res == 0)
1260         {
1261             cp_error();
1262         }
1263         if (magic != CPT_MAGIC2)
1264         {
1265             return -1;
1266         }
1267     }
1268
1269     return 0;
1270 }
1271
1272 static int do_cpt_state(XDR* xd, int fflags, t_state* state, FILE* list)
1273 {
1274     int             ret    = 0;
1275     const StatePart part   = StatePart::microState;
1276     const int       sflags = state->flags;
1277     // If reading, state->natoms was probably just read, so
1278     // allocations need to be managed. If writing, this won't change
1279     // anything that matters.
1280     state_change_natoms(state, state->natoms);
1281     for (int i = 0; (i < estNR && ret == 0); i++)
1282     {
1283         if (fflags & (1 << i))
1284         {
1285             switch (i)
1286             {
1287                 case estLAMBDA:
1288                     ret = doRealArrayRef(
1289                             xd,
1290                             part,
1291                             i,
1292                             sflags,
1293                             gmx::arrayRefFromArray<real>(state->lambda.data(), state->lambda.size()),
1294                             list);
1295                     break;
1296                 case estFEPSTATE:
1297                     ret = do_cpte_int(xd, part, i, sflags, &state->fep_state, list);
1298                     break;
1299                 case estBOX: ret = do_cpte_matrix(xd, part, i, sflags, state->box, list); break;
1300                 case estBOX_REL:
1301                     ret = do_cpte_matrix(xd, part, i, sflags, state->box_rel, list);
1302                     break;
1303                 case estBOXV: ret = do_cpte_matrix(xd, part, i, sflags, state->boxv, list); break;
1304                 case estPRES_PREV:
1305                     ret = do_cpte_matrix(xd, part, i, sflags, state->pres_prev, list);
1306                     break;
1307                 case estSVIR_PREV:
1308                     ret = do_cpte_matrix(xd, part, i, sflags, state->svir_prev, list);
1309                     break;
1310                 case estFVIR_PREV:
1311                     ret = do_cpte_matrix(xd, part, i, sflags, state->fvir_prev, list);
1312                     break;
1313                 case estNH_XI:
1314                     ret = doVector<double>(xd, part, i, sflags, &state->nosehoover_xi, list);
1315                     break;
1316                 case estNH_VXI:
1317                     ret = doVector<double>(xd, part, i, sflags, &state->nosehoover_vxi, list);
1318                     break;
1319                 case estNHPRES_XI:
1320                     ret = doVector<double>(xd, part, i, sflags, &state->nhpres_xi, list);
1321                     break;
1322                 case estNHPRES_VXI:
1323                     ret = doVector<double>(xd, part, i, sflags, &state->nhpres_vxi, list);
1324                     break;
1325                 case estTHERM_INT:
1326                     ret = doVector<double>(xd, part, i, sflags, &state->therm_integral, list);
1327                     break;
1328                 case estBAROS_INT:
1329                     ret = do_cpte_double(xd, part, i, sflags, &state->baros_integral, list);
1330                     break;
1331                 case estVETA: ret = do_cpte_real(xd, part, i, sflags, &state->veta, list); break;
1332                 case estVOL0: ret = do_cpte_real(xd, part, i, sflags, &state->vol0, list); break;
1333                 case estX:
1334                     ret = doRvecVector(xd, part, i, sflags, &state->x, state->natoms, list);
1335                     break;
1336                 case estV:
1337                     ret = doRvecVector(xd, part, i, sflags, &state->v, state->natoms, list);
1338                     break;
1339                 /* The RNG entries are no longer written,
1340                  * the next 4 lines are only for reading old files.
1341                  * It's OK that three case statements fall through.
1342                  */
1343                 case estLD_RNG_NOTSUPPORTED:
1344                 case estLD_RNGI_NOTSUPPORTED:
1345                 case estMC_RNG_NOTSUPPORTED:
1346                 case estMC_RNGI_NOTSUPPORTED:
1347                     ret = do_cpte_ints(xd, part, i, sflags, 0, nullptr, list);
1348                     break;
1349                 case estDISRE_INITF:
1350                     ret = do_cpte_real(xd, part, i, sflags, &state->hist.disre_initf, list);
1351                     break;
1352                 case estDISRE_RM3TAV:
1353                     ret = do_cpte_n_reals(
1354                             xd, part, i, sflags, &state->hist.ndisrepairs, &state->hist.disre_rm3tav, list);
1355                     break;
1356                 case estORIRE_INITF:
1357                     ret = do_cpte_real(xd, part, i, sflags, &state->hist.orire_initf, list);
1358                     break;
1359                 case estORIRE_DTAV:
1360                     ret = do_cpte_n_reals(
1361                             xd, part, i, sflags, &state->hist.norire_Dtav, &state->hist.orire_Dtav, list);
1362                     break;
1363                 case estPULLCOMPREVSTEP:
1364                     ret = doVector<double>(xd, part, i, sflags, &state->pull_com_prev_step, list);
1365                     break;
1366                 default:
1367                     gmx_fatal(FARGS,
1368                               "Unknown state entry %d\n"
1369                               "You are reading a checkpoint file written by different code, which "
1370                               "is not supported",
1371                               i);
1372             }
1373         }
1374     }
1375     return ret;
1376 }
1377
1378 static int do_cpt_ekinstate(XDR* xd, int fflags, ekinstate_t* ekins, FILE* list)
1379 {
1380     int ret = 0;
1381
1382     const StatePart part = StatePart::kineticEnergy;
1383     for (int i = 0; (i < eeksNR && ret == 0); i++)
1384     {
1385         if (fflags & (1 << i))
1386         {
1387             switch (i)
1388             {
1389
1390                 case eeksEKIN_N:
1391                     ret = do_cpte_int(xd, part, i, fflags, &ekins->ekin_n, list);
1392                     break;
1393                 case eeksEKINH:
1394                     ret = do_cpte_matrices(xd, part, i, fflags, ekins->ekin_n, &ekins->ekinh, list);
1395                     break;
1396                 case eeksEKINF:
1397                     ret = do_cpte_matrices(xd, part, i, fflags, ekins->ekin_n, &ekins->ekinf, list);
1398                     break;
1399                 case eeksEKINO:
1400                     ret = do_cpte_matrices(xd, part, i, fflags, ekins->ekin_n, &ekins->ekinh_old, list);
1401                     break;
1402                 case eeksEKINTOTAL:
1403                     ret = do_cpte_matrix(xd, part, i, fflags, ekins->ekin_total, list);
1404                     break;
1405                 case eeksEKINSCALEF:
1406                     ret = doVector<double>(xd, part, i, fflags, &ekins->ekinscalef_nhc, list);
1407                     break;
1408                 case eeksVSCALE:
1409                     ret = doVector<double>(xd, part, i, fflags, &ekins->vscale_nhc, list);
1410                     break;
1411                 case eeksEKINSCALEH:
1412                     ret = doVector<double>(xd, part, i, fflags, &ekins->ekinscaleh_nhc, list);
1413                     break;
1414                 case eeksDEKINDL:
1415                     ret = do_cpte_real(xd, part, i, fflags, &ekins->dekindl, list);
1416                     break;
1417                 case eeksMVCOS: ret = do_cpte_real(xd, part, i, fflags, &ekins->mvcos, list); break;
1418                 default:
1419                     gmx_fatal(FARGS,
1420                               "Unknown ekin data state entry %d\n"
1421                               "You are probably reading a new checkpoint file with old code",
1422                               i);
1423             }
1424         }
1425     }
1426
1427     return ret;
1428 }
1429
1430
1431 static int do_cpt_swapstate(XDR* xd, gmx_bool bRead, int eSwapCoords, swaphistory_t* swapstate, FILE* list)
1432 {
1433     int swap_cpt_version = 2;
1434
1435     if (eSwapCoords == eswapNO)
1436     {
1437         return 0;
1438     }
1439
1440     swapstate->bFromCpt    = bRead;
1441     swapstate->eSwapCoords = eSwapCoords;
1442
1443     do_cpt_int_err(xd, "swap checkpoint version", &swap_cpt_version, list);
1444     if (bRead && swap_cpt_version < 2)
1445     {
1446         gmx_fatal(FARGS,
1447                   "Cannot read checkpoint files that were written with old versions"
1448                   "of the ion/water position swapping protocol.\n");
1449     }
1450
1451     do_cpt_int_err(xd, "swap coupling steps", &swapstate->nAverage, list);
1452
1453     /* When reading, init_swapcoords has not been called yet,
1454      * so we have to allocate memory first. */
1455     do_cpt_int_err(xd, "number of ion types", &swapstate->nIonTypes, list);
1456     if (bRead)
1457     {
1458         snew(swapstate->ionType, swapstate->nIonTypes);
1459     }
1460
1461     for (int ic = 0; ic < eCompNR; ic++)
1462     {
1463         for (int ii = 0; ii < swapstate->nIonTypes; ii++)
1464         {
1465             swapstateIons_t* gs = &swapstate->ionType[ii];
1466
1467             if (bRead)
1468             {
1469                 do_cpt_int_err(xd, "swap requested atoms", &gs->nMolReq[ic], list);
1470             }
1471             else
1472             {
1473                 do_cpt_int_err(xd, "swap requested atoms p", gs->nMolReq_p[ic], list);
1474             }
1475
1476             if (bRead)
1477             {
1478                 do_cpt_int_err(xd, "swap influx net", &gs->inflow_net[ic], list);
1479             }
1480             else
1481             {
1482                 do_cpt_int_err(xd, "swap influx net p", gs->inflow_net_p[ic], list);
1483             }
1484
1485             if (bRead && (nullptr == gs->nMolPast[ic]))
1486             {
1487                 snew(gs->nMolPast[ic], swapstate->nAverage);
1488             }
1489
1490             for (int j = 0; j < swapstate->nAverage; j++)
1491             {
1492                 if (bRead)
1493                 {
1494                     do_cpt_int_err(xd, "swap past atom counts", &gs->nMolPast[ic][j], list);
1495                 }
1496                 else
1497                 {
1498                     do_cpt_int_err(xd, "swap past atom counts p", &gs->nMolPast_p[ic][j], list);
1499                 }
1500             }
1501         }
1502     }
1503
1504     /* Ion flux per channel */
1505     for (int ic = 0; ic < eChanNR; ic++)
1506     {
1507         for (int ii = 0; ii < swapstate->nIonTypes; ii++)
1508         {
1509             swapstateIons_t* gs = &swapstate->ionType[ii];
1510
1511             if (bRead)
1512             {
1513                 do_cpt_int_err(xd, "channel flux A->B", &gs->fluxfromAtoB[ic], list);
1514             }
1515             else
1516             {
1517                 do_cpt_int_err(xd, "channel flux A->B p", gs->fluxfromAtoB_p[ic], list);
1518             }
1519         }
1520     }
1521
1522     /* Ion flux leakage */
1523     if (bRead)
1524     {
1525         do_cpt_int_err(xd, "flux leakage", &swapstate->fluxleak, list);
1526     }
1527     else
1528     {
1529         do_cpt_int_err(xd, "flux leakage", swapstate->fluxleak_p, list);
1530     }
1531
1532     /* Ion history */
1533     for (int ii = 0; ii < swapstate->nIonTypes; ii++)
1534     {
1535         swapstateIons_t* gs = &swapstate->ionType[ii];
1536
1537         do_cpt_int_err(xd, "number of ions", &gs->nMol, list);
1538
1539         if (bRead)
1540         {
1541             snew(gs->channel_label, gs->nMol);
1542             snew(gs->comp_from, gs->nMol);
1543         }
1544
1545         do_cpt_u_chars(xd, "channel history", gs->nMol, gs->channel_label, list);
1546         do_cpt_u_chars(xd, "domain history", gs->nMol, gs->comp_from, list);
1547     }
1548
1549     /* Save the last known whole positions to checkpoint
1550      * file to be able to also make multimeric channels whole in PBC */
1551     do_cpt_int_err(xd, "Ch0 atoms", &swapstate->nat[eChan0], list);
1552     do_cpt_int_err(xd, "Ch1 atoms", &swapstate->nat[eChan1], list);
1553     if (bRead)
1554     {
1555         snew(swapstate->xc_old_whole[eChan0], swapstate->nat[eChan0]);
1556         snew(swapstate->xc_old_whole[eChan1], swapstate->nat[eChan1]);
1557         do_cpt_n_rvecs_err(xd, "Ch0 whole x", swapstate->nat[eChan0], swapstate->xc_old_whole[eChan0], list);
1558         do_cpt_n_rvecs_err(xd, "Ch1 whole x", swapstate->nat[eChan1], swapstate->xc_old_whole[eChan1], list);
1559     }
1560     else
1561     {
1562         do_cpt_n_rvecs_err(
1563                 xd, "Ch0 whole x", swapstate->nat[eChan0], *swapstate->xc_old_whole_p[eChan0], list);
1564         do_cpt_n_rvecs_err(
1565                 xd, "Ch1 whole x", swapstate->nat[eChan1], *swapstate->xc_old_whole_p[eChan1], list);
1566     }
1567
1568     return 0;
1569 }
1570
1571
1572 static int do_cpt_enerhist(XDR* xd, gmx_bool bRead, int fflags, energyhistory_t* enerhist, FILE* list)
1573 {
1574     int ret = 0;
1575
1576     if (fflags == 0)
1577     {
1578         return ret;
1579     }
1580
1581     GMX_RELEASE_ASSERT(enerhist != nullptr, "With energy history, we need a valid enerhist pointer");
1582
1583     /* This is stored/read for backward compatibility */
1584     int energyHistoryNumEnergies = 0;
1585     if (bRead)
1586     {
1587         enerhist->nsteps     = 0;
1588         enerhist->nsum       = 0;
1589         enerhist->nsteps_sim = 0;
1590         enerhist->nsum_sim   = 0;
1591     }
1592     else if (enerhist != nullptr)
1593     {
1594         energyHistoryNumEnergies = enerhist->ener_sum_sim.size();
1595     }
1596
1597     delta_h_history_t* deltaH = enerhist->deltaHForeignLambdas.get();
1598     const StatePart    part   = StatePart::energyHistory;
1599     for (int i = 0; (i < eenhNR && ret == 0); i++)
1600     {
1601         if (fflags & (1 << i))
1602         {
1603             switch (i)
1604             {
1605                 case eenhENERGY_N:
1606                     ret = do_cpte_int(xd, part, i, fflags, &energyHistoryNumEnergies, list);
1607                     break;
1608                 case eenhENERGY_AVER:
1609                     ret = doVector<double>(xd, part, i, fflags, &enerhist->ener_ave, list);
1610                     break;
1611                 case eenhENERGY_SUM:
1612                     ret = doVector<double>(xd, part, i, fflags, &enerhist->ener_sum, list);
1613                     break;
1614                 case eenhENERGY_NSUM:
1615                     do_cpt_step_err(xd, eenh_names[i], &enerhist->nsum, list);
1616                     break;
1617                 case eenhENERGY_SUM_SIM:
1618                     ret = doVector<double>(xd, part, i, fflags, &enerhist->ener_sum_sim, list);
1619                     break;
1620                 case eenhENERGY_NSUM_SIM:
1621                     do_cpt_step_err(xd, eenh_names[i], &enerhist->nsum_sim, list);
1622                     break;
1623                 case eenhENERGY_NSTEPS:
1624                     do_cpt_step_err(xd, eenh_names[i], &enerhist->nsteps, list);
1625                     break;
1626                 case eenhENERGY_NSTEPS_SIM:
1627                     do_cpt_step_err(xd, eenh_names[i], &enerhist->nsteps_sim, list);
1628                     break;
1629                 case eenhENERGY_DELTA_H_NN:
1630                 {
1631                     int numDeltaH = 0;
1632                     if (!bRead && deltaH != nullptr)
1633                     {
1634                         numDeltaH = deltaH->dh.size();
1635                     }
1636                     do_cpt_int_err(xd, eenh_names[i], &numDeltaH, list);
1637                     if (bRead)
1638                     {
1639                         if (deltaH == nullptr)
1640                         {
1641                             enerhist->deltaHForeignLambdas = std::make_unique<delta_h_history_t>();
1642                             deltaH                         = enerhist->deltaHForeignLambdas.get();
1643                         }
1644                         deltaH->dh.resize(numDeltaH);
1645                         deltaH->start_lambda_set = FALSE;
1646                     }
1647                     break;
1648                 }
1649                 case eenhENERGY_DELTA_H_LIST:
1650                     for (auto dh : deltaH->dh)
1651                     {
1652                         ret = doVector<real>(xd, part, i, fflags, &dh, list);
1653                     }
1654                     break;
1655                 case eenhENERGY_DELTA_H_STARTTIME:
1656                     ret = do_cpte_double(xd, part, i, fflags, &(deltaH->start_time), list);
1657                     break;
1658                 case eenhENERGY_DELTA_H_STARTLAMBDA:
1659                     ret = do_cpte_double(xd, part, i, fflags, &(deltaH->start_lambda), list);
1660                     break;
1661                 default:
1662                     gmx_fatal(FARGS,
1663                               "Unknown energy history entry %d\n"
1664                               "You are probably reading a new checkpoint file with old code",
1665                               i);
1666             }
1667         }
1668     }
1669
1670     if ((fflags & (1 << eenhENERGY_SUM)) && !(fflags & (1 << eenhENERGY_SUM_SIM)))
1671     {
1672         /* Assume we have an old file format and copy sum to sum_sim */
1673         enerhist->ener_sum_sim = enerhist->ener_sum;
1674     }
1675
1676     if ((fflags & (1 << eenhENERGY_NSUM)) && !(fflags & (1 << eenhENERGY_NSTEPS)))
1677     {
1678         /* Assume we have an old file format and copy nsum to nsteps */
1679         enerhist->nsteps = enerhist->nsum;
1680     }
1681     if ((fflags & (1 << eenhENERGY_NSUM_SIM)) && !(fflags & (1 << eenhENERGY_NSTEPS_SIM)))
1682     {
1683         /* Assume we have an old file format and copy nsum to nsteps */
1684         enerhist->nsteps_sim = enerhist->nsum_sim;
1685     }
1686
1687     return ret;
1688 }
1689
1690 static int doCptPullCoordHist(XDR* xd, PullCoordinateHistory* pullCoordHist, const StatePart part, FILE* list)
1691 {
1692     int ret   = 0;
1693     int flags = 0;
1694
1695     flags |= ((1 << epullcoordh_VALUE_REF_SUM) | (1 << epullcoordh_VALUE_SUM)
1696               | (1 << epullcoordh_DR01_SUM) | (1 << epullcoordh_DR23_SUM)
1697               | (1 << epullcoordh_DR45_SUM) | (1 << epullcoordh_FSCAL_SUM));
1698
1699     for (int i = 0; i < epullcoordh_NR && ret == 0; i++)
1700     {
1701         switch (i)
1702         {
1703             case epullcoordh_VALUE_REF_SUM:
1704                 ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->valueRef), list);
1705                 break;
1706             case epullcoordh_VALUE_SUM:
1707                 ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->value), list);
1708                 break;
1709             case epullcoordh_DR01_SUM:
1710                 for (int j = 0; j < DIM && ret == 0; j++)
1711                 {
1712                     ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dr01[j]), list);
1713                 }
1714                 break;
1715             case epullcoordh_DR23_SUM:
1716                 for (int j = 0; j < DIM && ret == 0; j++)
1717                 {
1718                     ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dr23[j]), list);
1719                 }
1720                 break;
1721             case epullcoordh_DR45_SUM:
1722                 for (int j = 0; j < DIM && ret == 0; j++)
1723                 {
1724                     ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dr45[j]), list);
1725                 }
1726                 break;
1727             case epullcoordh_FSCAL_SUM:
1728                 ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->scalarForce), list);
1729                 break;
1730             case epullcoordh_DYNAX_SUM:
1731                 for (int j = 0; j < DIM && ret == 0; j++)
1732                 {
1733                     ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dynaX[j]), list);
1734                 }
1735                 break;
1736         }
1737     }
1738
1739     return ret;
1740 }
1741
1742 static int doCptPullGroupHist(XDR* xd, PullGroupHistory* pullGroupHist, const StatePart part, FILE* list)
1743 {
1744     int ret   = 0;
1745     int flags = 0;
1746
1747     flags |= ((1 << epullgrouph_X_SUM));
1748
1749     for (int i = 0; i < epullgrouph_NR; i++)
1750     {
1751         switch (i)
1752         {
1753             case epullgrouph_X_SUM:
1754                 for (int j = 0; j < DIM && ret == 0; j++)
1755                 {
1756                     ret = do_cpte_double(xd, part, i, flags, &(pullGroupHist->x[j]), list);
1757                 }
1758                 break;
1759         }
1760     }
1761
1762     return ret;
1763 }
1764
1765
1766 static int doCptPullHist(XDR* xd, gmx_bool bRead, int fflags, PullHistory* pullHist, const StatePart part, FILE* list)
1767 {
1768     int ret                       = 0;
1769     int pullHistoryNumCoordinates = 0;
1770     int pullHistoryNumGroups      = 0;
1771
1772     /* Retain the number of terms in the sum and the number of coordinates (used for writing
1773      * average pull forces and coordinates) in the pull history, in temporary variables,
1774      * in case they cannot be read from the checkpoint, in order to have backward compatibility */
1775     if (bRead)
1776     {
1777         pullHist->numValuesInXSum = 0;
1778         pullHist->numValuesInFSum = 0;
1779     }
1780     else if (pullHist != nullptr)
1781     {
1782         pullHistoryNumCoordinates = pullHist->pullCoordinateSums.size();
1783         pullHistoryNumGroups      = pullHist->pullGroupSums.size();
1784     }
1785     else
1786     {
1787         GMX_RELEASE_ASSERT(fflags == 0, "Without pull history, all flags should be off");
1788     }
1789
1790     for (int i = 0; (i < epullhNR && ret == 0); i++)
1791     {
1792         if (fflags & (1 << i))
1793         {
1794             switch (i)
1795             {
1796                 case epullhPULL_NUMCOORDINATES:
1797                     ret = do_cpte_int(xd, part, i, fflags, &pullHistoryNumCoordinates, list);
1798                     break;
1799                 case epullhPULL_NUMGROUPS:
1800                     do_cpt_int_err(xd, eenh_names[i], &pullHistoryNumGroups, list);
1801                     break;
1802                 case epullhPULL_NUMVALUESINXSUM:
1803                     do_cpt_int_err(xd, eenh_names[i], &pullHist->numValuesInXSum, list);
1804                     break;
1805                 case epullhPULL_NUMVALUESINFSUM:
1806                     do_cpt_int_err(xd, eenh_names[i], &pullHist->numValuesInFSum, list);
1807                     break;
1808                 default:
1809                     gmx_fatal(FARGS,
1810                               "Unknown pull history entry %d\n"
1811                               "You are probably reading a new checkpoint file with old code",
1812                               i);
1813             }
1814         }
1815     }
1816     if (bRead)
1817     {
1818         pullHist->pullCoordinateSums.resize(pullHistoryNumCoordinates);
1819         pullHist->pullGroupSums.resize(pullHistoryNumGroups);
1820     }
1821     if (pullHist->numValuesInXSum > 0 || pullHist->numValuesInFSum > 0)
1822     {
1823         for (size_t i = 0; i < pullHist->pullCoordinateSums.size() && ret == 0; i++)
1824         {
1825             ret = doCptPullCoordHist(xd, &(pullHist->pullCoordinateSums[i]), part, list);
1826         }
1827         for (size_t i = 0; i < pullHist->pullGroupSums.size() && ret == 0; i++)
1828         {
1829             ret = doCptPullGroupHist(xd, &(pullHist->pullGroupSums[i]), part, list);
1830         }
1831     }
1832
1833     return ret;
1834 }
1835
1836 static int do_cpt_df_hist(XDR* xd, int fflags, int nlambda, df_history_t** dfhistPtr, FILE* list)
1837 {
1838     int ret = 0;
1839
1840     if (fflags == 0)
1841     {
1842         return 0;
1843     }
1844
1845     if (*dfhistPtr == nullptr)
1846     {
1847         snew(*dfhistPtr, 1);
1848         (*dfhistPtr)->nlambda = nlambda;
1849         init_df_history(*dfhistPtr, nlambda);
1850     }
1851     df_history_t* dfhist = *dfhistPtr;
1852
1853     const StatePart part = StatePart::freeEnergyHistory;
1854     for (int i = 0; (i < edfhNR && ret == 0); i++)
1855     {
1856         if (fflags & (1 << i))
1857         {
1858             switch (i)
1859             {
1860                 case edfhBEQUIL:
1861                     ret = do_cpte_bool(xd, part, i, fflags, &dfhist->bEquil, list);
1862                     break;
1863                 case edfhNATLAMBDA:
1864                     ret = do_cpte_ints(xd, part, i, fflags, nlambda, &dfhist->n_at_lam, list);
1865                     break;
1866                 case edfhWLHISTO:
1867                     ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->wl_histo, list);
1868                     break;
1869                 case edfhWLDELTA:
1870                     ret = do_cpte_real(xd, part, i, fflags, &dfhist->wl_delta, list);
1871                     break;
1872                 case edfhSUMWEIGHTS:
1873                     ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->sum_weights, list);
1874                     break;
1875                 case edfhSUMDG:
1876                     ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->sum_dg, list);
1877                     break;
1878                 case edfhSUMMINVAR:
1879                     ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->sum_minvar, list);
1880                     break;
1881                 case edfhSUMVAR:
1882                     ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->sum_variance, list);
1883                     break;
1884                 case edfhACCUMP:
1885                     ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->accum_p, list);
1886                     break;
1887                 case edfhACCUMM:
1888                     ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->accum_m, list);
1889                     break;
1890                 case edfhACCUMP2:
1891                     ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->accum_p2, list);
1892                     break;
1893                 case edfhACCUMM2:
1894                     ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->accum_m2, list);
1895                     break;
1896                 case edfhTIJ:
1897                     ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->Tij, list);
1898                     break;
1899                 case edfhTIJEMP:
1900                     ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->Tij_empirical, list);
1901                     break;
1902
1903                 default:
1904                     gmx_fatal(FARGS,
1905                               "Unknown df history entry %d\n"
1906                               "You are probably reading a new checkpoint file with old code",
1907                               i);
1908             }
1909         }
1910     }
1911
1912     return ret;
1913 }
1914
1915
1916 /* This function stores the last whole configuration of the reference and
1917  * average structure in the .cpt file
1918  */
1919 static int do_cpt_EDstate(XDR* xd, gmx_bool bRead, int nED, edsamhistory_t* EDstate, FILE* list)
1920 {
1921     if (nED == 0)
1922     {
1923         return 0;
1924     }
1925
1926     EDstate->bFromCpt = bRead;
1927     EDstate->nED      = nED;
1928
1929     /* When reading, init_edsam has not been called yet,
1930      * so we have to allocate memory first. */
1931     if (bRead)
1932     {
1933         snew(EDstate->nref, EDstate->nED);
1934         snew(EDstate->old_sref, EDstate->nED);
1935         snew(EDstate->nav, EDstate->nED);
1936         snew(EDstate->old_sav, EDstate->nED);
1937     }
1938
1939     /* Read/write the last whole conformation of SREF and SAV for each ED dataset (usually only one) */
1940     for (int i = 0; i < EDstate->nED; i++)
1941     {
1942         char buf[STRLEN];
1943
1944         /* Reference structure SREF */
1945         sprintf(buf, "ED%d # of atoms in reference structure", i + 1);
1946         do_cpt_int_err(xd, buf, &EDstate->nref[i], list);
1947         sprintf(buf, "ED%d x_ref", i + 1);
1948         if (bRead)
1949         {
1950             snew(EDstate->old_sref[i], EDstate->nref[i]);
1951             do_cpt_n_rvecs_err(xd, buf, EDstate->nref[i], EDstate->old_sref[i], list);
1952         }
1953         else
1954         {
1955             do_cpt_n_rvecs_err(xd, buf, EDstate->nref[i], EDstate->old_sref_p[i], list);
1956         }
1957
1958         /* Average structure SAV */
1959         sprintf(buf, "ED%d # of atoms in average structure", i + 1);
1960         do_cpt_int_err(xd, buf, &EDstate->nav[i], list);
1961         sprintf(buf, "ED%d x_av", i + 1);
1962         if (bRead)
1963         {
1964             snew(EDstate->old_sav[i], EDstate->nav[i]);
1965             do_cpt_n_rvecs_err(xd, buf, EDstate->nav[i], EDstate->old_sav[i], list);
1966         }
1967         else
1968         {
1969             do_cpt_n_rvecs_err(xd, buf, EDstate->nav[i], EDstate->old_sav_p[i], list);
1970         }
1971     }
1972
1973     return 0;
1974 }
1975
1976 static int do_cpt_correlation_grid(XDR*                         xd,
1977                                    gmx_bool                     bRead,
1978                                    gmx_unused int               fflags,
1979                                    gmx::CorrelationGridHistory* corrGrid,
1980                                    FILE*                        list,
1981                                    int                          eawhh)
1982 {
1983     int ret = 0;
1984
1985     do_cpt_int_err(xd, eawhh_names[eawhh], &(corrGrid->numCorrelationTensors), list);
1986     do_cpt_int_err(xd, eawhh_names[eawhh], &(corrGrid->tensorSize), list);
1987     do_cpt_int_err(xd, eawhh_names[eawhh], &(corrGrid->blockDataListSize), list);
1988
1989     if (bRead)
1990     {
1991         initCorrelationGridHistory(
1992                 corrGrid, corrGrid->numCorrelationTensors, corrGrid->tensorSize, corrGrid->blockDataListSize);
1993     }
1994
1995     for (gmx::CorrelationBlockDataHistory& blockData : corrGrid->blockDataBuffer)
1996     {
1997         do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumWeight), list);
1998         do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumSquareWeight), list);
1999         do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumWeightX), list);
2000         do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumWeightY), list);
2001         do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksSquareBlockWeight), list);
2002         do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksBlockSquareWeight), list);
2003         do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksBlockWeightBlockWeightX), list);
2004         do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksBlockWeightBlockWeightY), list);
2005         do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockLength), list);
2006         do_cpt_int_err(xd, eawhh_names[eawhh], &(blockData.previousBlockIndex), list);
2007         do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.correlationIntegral), list);
2008     }
2009
2010     return ret;
2011 }
2012
2013 static int do_cpt_awh_bias(XDR* xd, gmx_bool bRead, int fflags, gmx::AwhBiasHistory* biasHistory, FILE* list)
2014 {
2015     int ret = 0;
2016
2017     gmx::AwhBiasStateHistory* state = &biasHistory->state;
2018     for (int i = 0; (i < eawhhNR && ret == 0); i++)
2019     {
2020         if (fflags & (1 << i))
2021         {
2022             switch (i)
2023             {
2024                 case eawhhIN_INITIAL:
2025                     do_cpt_bool_err(xd, eawhh_names[i], &state->in_initial, list);
2026                     break;
2027                 case eawhhEQUILIBRATEHISTOGRAM:
2028                     do_cpt_bool_err(xd, eawhh_names[i], &state->equilibrateHistogram, list);
2029                     break;
2030                 case eawhhHISTSIZE:
2031                     do_cpt_double_err(xd, eawhh_names[i], &state->histSize, list);
2032                     break;
2033                 case eawhhNPOINTS:
2034                 {
2035                     int numPoints;
2036                     if (!bRead)
2037                     {
2038                         numPoints = biasHistory->pointState.size();
2039                     }
2040                     do_cpt_int_err(xd, eawhh_names[i], &numPoints, list);
2041                     if (bRead)
2042                     {
2043                         biasHistory->pointState.resize(numPoints);
2044                     }
2045                 }
2046                 break;
2047                 case eawhhCOORDPOINT:
2048                     for (auto& psh : biasHistory->pointState)
2049                     {
2050                         do_cpt_double_err(xd, eawhh_names[i], &psh.target, list);
2051                         do_cpt_double_err(xd, eawhh_names[i], &psh.free_energy, list);
2052                         do_cpt_double_err(xd, eawhh_names[i], &psh.bias, list);
2053                         do_cpt_double_err(xd, eawhh_names[i], &psh.weightsum_iteration, list);
2054                         do_cpt_double_err(xd, eawhh_names[i], &psh.weightsum_covering, list);
2055                         do_cpt_double_err(xd, eawhh_names[i], &psh.weightsum_tot, list);
2056                         do_cpt_double_err(xd, eawhh_names[i], &psh.weightsum_ref, list);
2057                         do_cpt_step_err(xd, eawhh_names[i], &psh.last_update_index, list);
2058                         do_cpt_double_err(xd, eawhh_names[i], &psh.log_pmfsum, list);
2059                         do_cpt_double_err(xd, eawhh_names[i], &psh.visits_iteration, list);
2060                         do_cpt_double_err(xd, eawhh_names[i], &psh.visits_tot, list);
2061                     }
2062                     break;
2063                 case eawhhUMBRELLAGRIDPOINT:
2064                     do_cpt_int_err(xd, eawhh_names[i], &(state->umbrellaGridpoint), list);
2065                     break;
2066                 case eawhhUPDATELIST:
2067                     do_cpt_int_err(xd, eawhh_names[i], &(state->origin_index_updatelist), list);
2068                     do_cpt_int_err(xd, eawhh_names[i], &(state->end_index_updatelist), list);
2069                     break;
2070                 case eawhhLOGSCALEDSAMPLEWEIGHT:
2071                     do_cpt_double_err(xd, eawhh_names[i], &(state->logScaledSampleWeight), list);
2072                     do_cpt_double_err(xd, eawhh_names[i], &(state->maxLogScaledSampleWeight), list);
2073                     break;
2074                 case eawhhNUMUPDATES:
2075                     do_cpt_step_err(xd, eawhh_names[i], &(state->numUpdates), list);
2076                     break;
2077                 case eawhhFORCECORRELATIONGRID:
2078                     ret = do_cpt_correlation_grid(
2079                             xd, bRead, fflags, &biasHistory->forceCorrelationGrid, list, i);
2080                     break;
2081                 default: gmx_fatal(FARGS, "Unknown awh history entry %d\n", i);
2082             }
2083         }
2084     }
2085
2086     return ret;
2087 }
2088
2089 static int do_cpt_awh(XDR* xd, gmx_bool bRead, int fflags, gmx::AwhHistory* awhHistory, FILE* list)
2090 {
2091     int ret = 0;
2092
2093     if (fflags != 0)
2094     {
2095         std::shared_ptr<gmx::AwhHistory> awhHistoryLocal;
2096
2097         if (awhHistory == nullptr)
2098         {
2099             GMX_RELEASE_ASSERT(bRead,
2100                                "do_cpt_awh should not be called for writing without an AwhHistory");
2101
2102             awhHistoryLocal = std::make_shared<gmx::AwhHistory>();
2103             awhHistory      = awhHistoryLocal.get();
2104         }
2105
2106         /* To be able to read and write the AWH data several parameters determining the layout of the AWH structs need to be known
2107            (nbias, npoints, etc.). The best thing (?) would be to have these passed to this function. When writing to a checkpoint
2108            these parameters are available in awh_history (after calling init_awh_history). When reading from a checkpoint though, there
2109            is no initialized awh_history (it is initialized and set in this function). The AWH parameters have not always been read
2110            at the time when this function is called for reading so I don't know how to pass them as input. Here, this is solved by
2111            when writing a checkpoint, also storing parameters needed for future reading of the checkpoint.
2112
2113            Another issue is that some variables that AWH checkpoints don't have a registered enum and string (e.g. nbias below).
2114            One difficulty is the multilevel structure of the data which would need to be represented somehow. */
2115
2116         int numBias;
2117         if (!bRead)
2118         {
2119             numBias = awhHistory->bias.size();
2120         }
2121         do_cpt_int_err(xd, "awh_nbias", &numBias, list);
2122
2123         if (bRead)
2124         {
2125             awhHistory->bias.resize(numBias);
2126         }
2127         for (auto& bias : awhHistory->bias)
2128         {
2129             ret = do_cpt_awh_bias(xd, bRead, fflags, &bias, list);
2130             if (ret)
2131             {
2132                 return ret;
2133             }
2134         }
2135         do_cpt_double_err(xd, "awh_potential_offset", &awhHistory->potentialOffset, list);
2136     }
2137     return ret;
2138 }
2139
2140 static void do_cpt_mdmodules(int                           fileVersion,
2141                              t_fileio*                     checkpointFileHandle,
2142                              const gmx::MdModulesNotifier& mdModulesNotifier,
2143                              FILE*                         outputFile)
2144 {
2145     if (fileVersion >= cptv_MdModules)
2146     {
2147         gmx::FileIOXdrSerializer serializer(checkpointFileHandle);
2148         gmx::KeyValueTreeObject  mdModuleCheckpointParameterTree =
2149                 gmx::deserializeKeyValueTree(&serializer);
2150         if (outputFile)
2151         {
2152             gmx::TextWriter textWriter(outputFile);
2153             gmx::dumpKeyValueTree(&textWriter, mdModuleCheckpointParameterTree);
2154         }
2155         gmx::MdModulesCheckpointReadingDataOnMaster mdModuleCheckpointReadingDataOnMaster = {
2156             mdModuleCheckpointParameterTree, fileVersion
2157         };
2158         mdModulesNotifier.checkpointingNotifications_.notify(mdModuleCheckpointReadingDataOnMaster);
2159     }
2160 }
2161
2162 static int do_cpt_files(XDR* xd, gmx_bool bRead, std::vector<gmx_file_position_t>* outputfiles, FILE* list, int file_version)
2163 {
2164     gmx_off_t                   offset;
2165     gmx_off_t                   mask = 0xFFFFFFFFL;
2166     int                         offset_high, offset_low;
2167     std::array<char, CPTSTRLEN> buf;
2168     GMX_RELEASE_ASSERT(outputfiles, "Must have valid outputfiles");
2169
2170     // Ensure that reading pre-allocates outputfiles, while writing
2171     // writes what is already there.
2172     int nfiles = outputfiles->size();
2173     if (do_cpt_int(xd, "number of output files", &nfiles, list) != 0)
2174     {
2175         return -1;
2176     }
2177     if (bRead)
2178     {
2179         outputfiles->resize(nfiles);
2180     }
2181
2182     for (auto& outputfile : *outputfiles)
2183     {
2184         /* 64-bit XDR numbers are not portable, so it is stored as separate high/low fractions */
2185         if (bRead)
2186         {
2187             do_cpt_string_err(xd, "output filename", buf, list);
2188             std::copy(std::begin(buf), std::end(buf), std::begin(outputfile.filename));
2189
2190             if (do_cpt_int(xd, "file_offset_high", &offset_high, list) != 0)
2191             {
2192                 return -1;
2193             }
2194             if (do_cpt_int(xd, "file_offset_low", &offset_low, list) != 0)
2195             {
2196                 return -1;
2197             }
2198             outputfile.offset = (static_cast<gmx_off_t>(offset_high) << 32)
2199                                 | (static_cast<gmx_off_t>(offset_low) & mask);
2200         }
2201         else
2202         {
2203             do_cpt_string_err(xd, "output filename", outputfile.filename, list);
2204             /* writing */
2205             offset = outputfile.offset;
2206             if (offset == -1)
2207             {
2208                 offset_low  = -1;
2209                 offset_high = -1;
2210             }
2211             else
2212             {
2213                 offset_low  = static_cast<int>(offset & mask);
2214                 offset_high = static_cast<int>((offset >> 32) & mask);
2215             }
2216             if (do_cpt_int(xd, "file_offset_high", &offset_high, list) != 0)
2217             {
2218                 return -1;
2219             }
2220             if (do_cpt_int(xd, "file_offset_low", &offset_low, list) != 0)
2221             {
2222                 return -1;
2223             }
2224         }
2225         if (file_version >= 8)
2226         {
2227             if (do_cpt_int(xd, "file_checksum_size", &outputfile.checksumSize, list) != 0)
2228             {
2229                 return -1;
2230             }
2231             if (do_cpt_u_chars(xd, "file_checksum", outputfile.checksum.size(), outputfile.checksum.data(), list)
2232                 != 0)
2233             {
2234                 return -1;
2235             }
2236         }
2237         else
2238         {
2239             outputfile.checksumSize = -1;
2240         }
2241     }
2242     return 0;
2243 }
2244
2245 void write_checkpoint_data(t_fileio*                         fp,
2246                            CheckpointHeaderContents          headerContents,
2247                            gmx_bool                          bExpanded,
2248                            int                               elamstats,
2249                            t_state*                          state,
2250                            ObservablesHistory*               observablesHistory,
2251                            const gmx::MdModulesNotifier&     mdModulesNotifier,
2252                            std::vector<gmx_file_position_t>* outputfiles,
2253                            gmx::WriteCheckpointDataHolder*   modularSimulatorCheckpointData)
2254 {
2255     headerContents.flags_eks = 0;
2256     if (state->ekinstate.bUpToDate)
2257     {
2258         headerContents.flags_eks = ((1 << eeksEKIN_N) | (1 << eeksEKINH) | (1 << eeksEKINF)
2259                                     | (1 << eeksEKINO) | (1 << eeksEKINSCALEF) | (1 << eeksEKINSCALEH)
2260                                     | (1 << eeksVSCALE) | (1 << eeksDEKINDL) | (1 << eeksMVCOS));
2261     }
2262     headerContents.isModularSimulatorCheckpoint = !modularSimulatorCheckpointData->empty();
2263
2264     energyhistory_t* enerhist = observablesHistory->energyHistory.get();
2265     headerContents.flags_enh  = 0;
2266     if (enerhist != nullptr && (enerhist->nsum > 0 || enerhist->nsum_sim > 0))
2267     {
2268         headerContents.flags_enh |=
2269                 (1 << eenhENERGY_N) | (1 << eenhENERGY_NSTEPS) | (1 << eenhENERGY_NSTEPS_SIM);
2270         if (enerhist->nsum > 0)
2271         {
2272             headerContents.flags_enh |=
2273                     ((1 << eenhENERGY_AVER) | (1 << eenhENERGY_SUM) | (1 << eenhENERGY_NSUM));
2274         }
2275         if (enerhist->nsum_sim > 0)
2276         {
2277             headerContents.flags_enh |= ((1 << eenhENERGY_SUM_SIM) | (1 << eenhENERGY_NSUM_SIM));
2278         }
2279         if (enerhist->deltaHForeignLambdas != nullptr)
2280         {
2281             headerContents.flags_enh |=
2282                     ((1 << eenhENERGY_DELTA_H_NN) | (1 << eenhENERGY_DELTA_H_LIST)
2283                      | (1 << eenhENERGY_DELTA_H_STARTTIME) | (1 << eenhENERGY_DELTA_H_STARTLAMBDA));
2284         }
2285     }
2286
2287     PullHistory* pullHist           = observablesHistory->pullHistory.get();
2288     headerContents.flagsPullHistory = 0;
2289     if (pullHist != nullptr && (pullHist->numValuesInXSum > 0 || pullHist->numValuesInFSum > 0))
2290     {
2291         headerContents.flagsPullHistory |= (1 << epullhPULL_NUMCOORDINATES);
2292         headerContents.flagsPullHistory |= ((1 << epullhPULL_NUMGROUPS) | (1 << epullhPULL_NUMVALUESINXSUM)
2293                                             | (1 << epullhPULL_NUMVALUESINFSUM));
2294     }
2295
2296     headerContents.flags_dfh = 0;
2297     if (bExpanded)
2298     {
2299         headerContents.flags_dfh = ((1 << edfhBEQUIL) | (1 << edfhNATLAMBDA) | (1 << edfhSUMWEIGHTS)
2300                                     | (1 << edfhSUMDG) | (1 << edfhTIJ) | (1 << edfhTIJEMP));
2301         if (EWL(elamstats))
2302         {
2303             headerContents.flags_dfh |= ((1 << edfhWLDELTA) | (1 << edfhWLHISTO));
2304         }
2305         if ((elamstats == elamstatsMINVAR) || (elamstats == elamstatsBARKER)
2306             || (elamstats == elamstatsMETROPOLIS))
2307         {
2308             headerContents.flags_dfh |= ((1 << edfhACCUMP) | (1 << edfhACCUMM) | (1 << edfhACCUMP2)
2309                                          | (1 << edfhACCUMM2) | (1 << edfhSUMMINVAR) | (1 << edfhSUMVAR));
2310         }
2311     }
2312
2313     headerContents.flags_awhh = 0;
2314     if (state->awhHistory != nullptr && !state->awhHistory->bias.empty())
2315     {
2316         headerContents.flags_awhh |=
2317                 ((1 << eawhhIN_INITIAL) | (1 << eawhhEQUILIBRATEHISTOGRAM) | (1 << eawhhHISTSIZE)
2318                  | (1 << eawhhNPOINTS) | (1 << eawhhCOORDPOINT) | (1 << eawhhUMBRELLAGRIDPOINT)
2319                  | (1 << eawhhUPDATELIST) | (1 << eawhhLOGSCALEDSAMPLEWEIGHT)
2320                  | (1 << eawhhNUMUPDATES) | (1 << eawhhFORCECORRELATIONGRID));
2321     }
2322
2323     do_cpt_header(gmx_fio_getxdr(fp), FALSE, nullptr, &headerContents);
2324
2325     if ((do_cpt_state(gmx_fio_getxdr(fp), state->flags, state, nullptr) < 0)
2326         || (do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents.flags_eks, &state->ekinstate, nullptr) < 0)
2327         || (do_cpt_enerhist(gmx_fio_getxdr(fp), FALSE, headerContents.flags_enh, enerhist, nullptr) < 0)
2328         || (doCptPullHist(gmx_fio_getxdr(fp), FALSE, headerContents.flagsPullHistory, pullHist, StatePart::pullHistory, nullptr)
2329             < 0)
2330         || (do_cpt_df_hist(gmx_fio_getxdr(fp), headerContents.flags_dfh, headerContents.nlambda, &state->dfhist, nullptr)
2331             < 0)
2332         || (do_cpt_EDstate(
2333                     gmx_fio_getxdr(fp), FALSE, headerContents.nED, observablesHistory->edsamHistory.get(), nullptr)
2334             < 0)
2335         || (do_cpt_awh(gmx_fio_getxdr(fp), FALSE, headerContents.flags_awhh, state->awhHistory.get(), nullptr) < 0)
2336         || (do_cpt_swapstate(gmx_fio_getxdr(fp),
2337                              FALSE,
2338                              headerContents.eSwapCoords,
2339                              observablesHistory->swapHistory.get(),
2340                              nullptr)
2341             < 0)
2342         || (do_cpt_files(gmx_fio_getxdr(fp), FALSE, outputfiles, nullptr, headerContents.file_version) < 0))
2343     {
2344         gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
2345     }
2346
2347     // Checkpointing MdModules
2348     {
2349         gmx::KeyValueTreeBuilder          builder;
2350         gmx::MdModulesWriteCheckpointData mdModulesWriteCheckpoint = { builder.rootObject(),
2351                                                                        headerContents.file_version };
2352         mdModulesNotifier.checkpointingNotifications_.notify(mdModulesWriteCheckpoint);
2353         auto                     tree = builder.build();
2354         gmx::FileIOXdrSerializer serializer(fp);
2355         gmx::serializeKeyValueTree(tree, &serializer);
2356     }
2357
2358     // Checkpointing modular simulator
2359     {
2360         gmx::FileIOXdrSerializer serializer(fp);
2361         modularSimulatorCheckpointData->serialize(&serializer);
2362     }
2363
2364     do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
2365 #if GMX_FAHCORE
2366     /* Always FAH checkpoint immediately after a Gromacs checkpoint.
2367      *
2368      * Note that it is critical that we save a FAH checkpoint directly
2369      * after writing a Gromacs checkpoint.  If the program dies, either
2370      * by the machine powering off suddenly or the process being,
2371      * killed, FAH can recover files that have only appended data by
2372      * truncating them to the last recorded length.  The Gromacs
2373      * checkpoint does not just append data, it is fully rewritten each
2374      * time so a crash between moving the new Gromacs checkpoint file in
2375      * to place and writing a FAH checkpoint is not recoverable.  Thus
2376      * the time between these operations must be kept as short a
2377      * possible.
2378      */
2379     fcCheckpoint();
2380 #endif
2381 }
2382
2383 static void check_int(FILE* fplog, const char* type, int p, int f, gmx_bool* mm)
2384 {
2385     bool foundMismatch = (p != f);
2386     if (!foundMismatch)
2387     {
2388         return;
2389     }
2390     *mm = TRUE;
2391     if (fplog)
2392     {
2393         fprintf(fplog, "  %s mismatch,\n", type);
2394         fprintf(fplog, "    current program: %d\n", p);
2395         fprintf(fplog, "    checkpoint file: %d\n", f);
2396         fprintf(fplog, "\n");
2397     }
2398 }
2399
2400 static void check_string(FILE* fplog, const char* type, const char* p, const char* f, gmx_bool* mm)
2401 {
2402     bool foundMismatch = (std::strcmp(p, f) != 0);
2403     if (!foundMismatch)
2404     {
2405         return;
2406     }
2407     *mm = TRUE;
2408     if (fplog)
2409     {
2410         fprintf(fplog, "  %s mismatch,\n", type);
2411         fprintf(fplog, "    current program: %s\n", p);
2412         fprintf(fplog, "    checkpoint file: %s\n", f);
2413         fprintf(fplog, "\n");
2414     }
2415 }
2416
2417 static void check_match(FILE*                           fplog,
2418                         const t_commrec*                cr,
2419                         const ivec                      dd_nc,
2420                         const CheckpointHeaderContents& headerContents,
2421                         gmx_bool                        reproducibilityRequested)
2422 {
2423     /* Note that this check_string on the version will also print a message
2424      * when only the minor version differs. But we only print a warning
2425      * message further down with reproducibilityRequested=TRUE.
2426      */
2427     gmx_bool versionDiffers = FALSE;
2428     check_string(fplog, "Version", gmx_version(), headerContents.version, &versionDiffers);
2429
2430     gmx_bool precisionDiffers = FALSE;
2431     check_int(fplog, "Double prec.", GMX_DOUBLE, headerContents.double_prec, &precisionDiffers);
2432     if (precisionDiffers)
2433     {
2434         const char msg_precision_difference[] =
2435                 "You are continuing a simulation with a different precision. Not matching\n"
2436                 "mixed/double precision will lead to precision or performance loss.\n";
2437         if (fplog)
2438         {
2439             fprintf(fplog, "%s\n", msg_precision_difference);
2440         }
2441     }
2442
2443     gmx_bool mm = (versionDiffers || precisionDiffers);
2444
2445     if (reproducibilityRequested)
2446     {
2447         check_string(
2448                 fplog, "Program name", gmx::getProgramContext().fullBinaryPath(), headerContents.fprog, &mm);
2449
2450         check_int(fplog, "#ranks", cr->nnodes, headerContents.nnodes, &mm);
2451     }
2452
2453     if (cr->sizeOfDefaultCommunicator > 1 && reproducibilityRequested)
2454     {
2455         // TODO: These checks are incorrect (see redmine #3309)
2456         check_int(fplog, "#PME-ranks", cr->npmenodes, headerContents.npme, &mm);
2457
2458         int npp = cr->sizeOfDefaultCommunicator;
2459         if (cr->npmenodes >= 0)
2460         {
2461             npp -= cr->npmenodes;
2462         }
2463         int npp_f = headerContents.nnodes;
2464         if (headerContents.npme >= 0)
2465         {
2466             npp_f -= headerContents.npme;
2467         }
2468         if (npp == npp_f)
2469         {
2470             check_int(fplog, "#DD-cells[x]", dd_nc[XX], headerContents.dd_nc[XX], &mm);
2471             check_int(fplog, "#DD-cells[y]", dd_nc[YY], headerContents.dd_nc[YY], &mm);
2472             check_int(fplog, "#DD-cells[z]", dd_nc[ZZ], headerContents.dd_nc[ZZ], &mm);
2473         }
2474     }
2475
2476     if (mm)
2477     {
2478         /* Gromacs should be able to continue from checkpoints between
2479          * different patch level versions, but we do not guarantee
2480          * compatibility between different major/minor versions - check this.
2481          */
2482         int gmx_major;
2483         int cpt_major;
2484         sscanf(gmx_version(), "%5d", &gmx_major);
2485         int      ret                 = sscanf(headerContents.version, "%5d", &cpt_major);
2486         gmx_bool majorVersionDiffers = (ret < 1 || gmx_major != cpt_major);
2487
2488         const char msg_major_version_difference[] =
2489                 "The current GROMACS major version is not identical to the one that\n"
2490                 "generated the checkpoint file. In principle GROMACS does not support\n"
2491                 "continuation from checkpoints between different versions, so we advise\n"
2492                 "against this. If you still want to try your luck we recommend that you use\n"
2493                 "the -noappend flag to keep your output files from the two versions separate.\n"
2494                 "This might also work around errors where the output fields in the energy\n"
2495                 "file have changed between the different versions.\n";
2496
2497         const char msg_mismatch_notice[] =
2498                 "GROMACS patchlevel, binary or parallel settings differ from previous run.\n"
2499                 "Continuation is exact, but not guaranteed to be binary identical.\n";
2500
2501         if (majorVersionDiffers)
2502         {
2503             if (fplog)
2504             {
2505                 fprintf(fplog, "%s\n", msg_major_version_difference);
2506             }
2507         }
2508         else if (reproducibilityRequested)
2509         {
2510             /* Major & minor versions match at least, but something is different. */
2511             if (fplog)
2512             {
2513                 fprintf(fplog, "%s\n", msg_mismatch_notice);
2514             }
2515         }
2516     }
2517 }
2518
2519 static void read_checkpoint(const char*                    fn,
2520                             t_fileio*                      logfio,
2521                             const t_commrec*               cr,
2522                             const ivec                     dd_nc,
2523                             int                            eIntegrator,
2524                             int*                           init_fep_state,
2525                             CheckpointHeaderContents*      headerContents,
2526                             t_state*                       state,
2527                             ObservablesHistory*            observablesHistory,
2528                             gmx_bool                       reproducibilityRequested,
2529                             const gmx::MdModulesNotifier&  mdModulesNotifier,
2530                             gmx::ReadCheckpointDataHolder* modularSimulatorCheckpointData,
2531                             bool                           useModularSimulator)
2532 {
2533     t_fileio* fp;
2534     char      buf[STEPSTRSIZE];
2535     int       ret;
2536
2537     fp = gmx_fio_open(fn, "r");
2538     do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, headerContents);
2539
2540     // If we are appending, then we don't want write to the open log
2541     // file because we still need to compute a checksum for it. In
2542     // that case, the filehandle will be nullptr. Otherwise, we report
2543     // to the new log file about the checkpoint file that we are
2544     // reading from.
2545     FILE* fplog = gmx_fio_getfp(logfio);
2546     if (fplog)
2547     {
2548         fprintf(fplog, "\n");
2549         fprintf(fplog, "Reading checkpoint file %s\n", fn);
2550         fprintf(fplog, "  file generated by:     %s\n", headerContents->fprog);
2551         fprintf(fplog, "  file generated at:     %s\n", headerContents->ftime);
2552         fprintf(fplog, "  GROMACS double prec.:  %d\n", headerContents->double_prec);
2553         fprintf(fplog, "  simulation part #:     %d\n", headerContents->simulation_part);
2554         fprintf(fplog, "  step:                  %s\n", gmx_step_str(headerContents->step, buf));
2555         fprintf(fplog, "  time:                  %f\n", headerContents->t);
2556         fprintf(fplog, "\n");
2557     }
2558
2559     if (headerContents->natoms != state->natoms)
2560     {
2561         gmx_fatal(FARGS,
2562                   "Checkpoint file is for a system of %d atoms, while the current system consists "
2563                   "of %d atoms",
2564                   headerContents->natoms,
2565                   state->natoms);
2566     }
2567     if (headerContents->ngtc != state->ngtc)
2568     {
2569         gmx_fatal(FARGS,
2570                   "Checkpoint file is for a system of %d T-coupling groups, while the current "
2571                   "system consists of %d T-coupling groups",
2572                   headerContents->ngtc,
2573                   state->ngtc);
2574     }
2575     if (headerContents->nnhpres != state->nnhpres)
2576     {
2577         gmx_fatal(FARGS,
2578                   "Checkpoint file is for a system of %d NH-pressure-coupling variables, while the "
2579                   "current system consists of %d NH-pressure-coupling variables",
2580                   headerContents->nnhpres,
2581                   state->nnhpres);
2582     }
2583
2584     int nlambdaHistory = (state->dfhist ? state->dfhist->nlambda : 0);
2585     if (headerContents->nlambda != nlambdaHistory)
2586     {
2587         gmx_fatal(FARGS,
2588                   "Checkpoint file is for a system with %d lambda states, while the current system "
2589                   "consists of %d lambda states",
2590                   headerContents->nlambda,
2591                   nlambdaHistory);
2592     }
2593
2594     init_gtc_state(state,
2595                    state->ngtc,
2596                    state->nnhpres,
2597                    headerContents->nhchainlength); /* need to keep this here to keep the tpr format working */
2598     /* write over whatever was read; we use the number of Nose-Hoover chains from the checkpoint */
2599
2600     if (headerContents->eIntegrator != eIntegrator)
2601     {
2602         gmx_fatal(FARGS,
2603                   "Cannot change integrator during a checkpoint restart. Perhaps you should make a "
2604                   "new .tpr with grompp -f new.mdp -t %s",
2605                   fn);
2606     }
2607
2608     // For modular simulator, no state object is populated, so we cannot do this check here!
2609     if (headerContents->flags_state != state->flags && !useModularSimulator)
2610     {
2611         gmx_fatal(FARGS,
2612                   "Cannot change a simulation algorithm during a checkpoint restart. Perhaps you "
2613                   "should make a new .tpr with grompp -f new.mdp -t %s",
2614                   fn);
2615     }
2616
2617     GMX_RELEASE_ASSERT(!(headerContents->isModularSimulatorCheckpoint && !useModularSimulator),
2618                        "Checkpoint file was written by modular simulator, but the current "
2619                        "simulation uses the legacy simulator.\n\n"
2620                        "Try the following steps:\n"
2621                        "1. Make sure the GMX_DISABLE_MODULAR_SIMULATOR environment variable is not "
2622                        "set to return to the default behavior. Retry running the simulation.\n"
2623                        "2. If the problem persists, set the environment variable "
2624                        "GMX_USE_MODULAR_SIMULATOR=ON to overwrite the default behavior and use "
2625                        "modular simulator for all implemented use cases.");
2626     GMX_RELEASE_ASSERT(!(!headerContents->isModularSimulatorCheckpoint && useModularSimulator),
2627                        "Checkpoint file was written by legacy simulator, but the current "
2628                        "simulation uses the modular simulator.\n\n"
2629                        "Try the following steps:\n"
2630                        "1. Make sure the GMX_USE_MODULAR_SIMULATOR environment variable is not set "
2631                        "to return to the default behavior. Retry running the simulation.\n"
2632                        "2. If the problem persists, set the environment variable "
2633                        "GMX_DISABLE_MODULAR_SIMULATOR=ON to overwrite the default behavior and use "
2634                        "legacy simulator for all implemented use cases.");
2635
2636     if (MASTER(cr))
2637     {
2638         check_match(fplog, cr, dd_nc, *headerContents, reproducibilityRequested);
2639     }
2640
2641     ret             = do_cpt_state(gmx_fio_getxdr(fp), headerContents->flags_state, state, nullptr);
2642     *init_fep_state = state->fep_state; /* there should be a better way to do this than setting it
2643                                            here. Investigate for 5.0. */
2644     if (ret)
2645     {
2646         cp_error();
2647     }
2648     ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents->flags_eks, &state->ekinstate, nullptr);
2649     if (ret)
2650     {
2651         cp_error();
2652     }
2653     state->ekinstate.hasReadEkinState = (((headerContents->flags_eks & (1 << eeksEKINH)) != 0)
2654                                          || ((headerContents->flags_eks & (1 << eeksEKINF)) != 0)
2655                                          || ((headerContents->flags_eks & (1 << eeksEKINO)) != 0)
2656                                          || (((headerContents->flags_eks & (1 << eeksEKINSCALEF))
2657                                               | (headerContents->flags_eks & (1 << eeksEKINSCALEH))
2658                                               | (headerContents->flags_eks & (1 << eeksVSCALE)))
2659                                              != 0));
2660
2661     if (headerContents->flags_enh && observablesHistory->energyHistory == nullptr)
2662     {
2663         observablesHistory->energyHistory = std::make_unique<energyhistory_t>();
2664     }
2665     ret = do_cpt_enerhist(
2666             gmx_fio_getxdr(fp), TRUE, headerContents->flags_enh, observablesHistory->energyHistory.get(), nullptr);
2667     if (ret)
2668     {
2669         cp_error();
2670     }
2671
2672     if (headerContents->flagsPullHistory)
2673     {
2674         if (observablesHistory->pullHistory == nullptr)
2675         {
2676             observablesHistory->pullHistory = std::make_unique<PullHistory>();
2677         }
2678         ret = doCptPullHist(gmx_fio_getxdr(fp),
2679                             TRUE,
2680                             headerContents->flagsPullHistory,
2681                             observablesHistory->pullHistory.get(),
2682                             StatePart::pullHistory,
2683                             nullptr);
2684         if (ret)
2685         {
2686             cp_error();
2687         }
2688     }
2689
2690     if (headerContents->file_version < 6)
2691     {
2692         gmx_fatal(FARGS,
2693                   "Continuing from checkpoint files written before GROMACS 4.5 is not supported");
2694     }
2695
2696     ret = do_cpt_df_hist(
2697             gmx_fio_getxdr(fp), headerContents->flags_dfh, headerContents->nlambda, &state->dfhist, nullptr);
2698     if (ret)
2699     {
2700         cp_error();
2701     }
2702
2703     if (headerContents->nED > 0 && observablesHistory->edsamHistory == nullptr)
2704     {
2705         observablesHistory->edsamHistory = std::make_unique<edsamhistory_t>(edsamhistory_t{});
2706     }
2707     ret = do_cpt_EDstate(
2708             gmx_fio_getxdr(fp), TRUE, headerContents->nED, observablesHistory->edsamHistory.get(), nullptr);
2709     if (ret)
2710     {
2711         cp_error();
2712     }
2713
2714     if (headerContents->flags_awhh != 0 && state->awhHistory == nullptr)
2715     {
2716         state->awhHistory = std::make_shared<gmx::AwhHistory>();
2717     }
2718     ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE, headerContents->flags_awhh, state->awhHistory.get(), nullptr);
2719     if (ret)
2720     {
2721         cp_error();
2722     }
2723
2724     if (headerContents->eSwapCoords != eswapNO && observablesHistory->swapHistory == nullptr)
2725     {
2726         observablesHistory->swapHistory = std::make_unique<swaphistory_t>(swaphistory_t{});
2727     }
2728     ret = do_cpt_swapstate(
2729             gmx_fio_getxdr(fp), TRUE, headerContents->eSwapCoords, observablesHistory->swapHistory.get(), nullptr);
2730     if (ret)
2731     {
2732         cp_error();
2733     }
2734
2735     std::vector<gmx_file_position_t> outputfiles;
2736     ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, nullptr, headerContents->file_version);
2737     if (ret)
2738     {
2739         cp_error();
2740     }
2741     do_cpt_mdmodules(headerContents->file_version, fp, mdModulesNotifier, nullptr);
2742     if (headerContents->file_version >= cptv_ModularSimulator)
2743     {
2744         gmx::FileIOXdrSerializer serializer(fp);
2745         modularSimulatorCheckpointData->deserialize(&serializer);
2746     }
2747     ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents->file_version);
2748     if (ret)
2749     {
2750         cp_error();
2751     }
2752     if (gmx_fio_close(fp) != 0)
2753     {
2754         gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
2755     }
2756 }
2757
2758
2759 void load_checkpoint(const char*                    fn,
2760                      t_fileio*                      logfio,
2761                      const t_commrec*               cr,
2762                      const ivec                     dd_nc,
2763                      t_inputrec*                    ir,
2764                      t_state*                       state,
2765                      ObservablesHistory*            observablesHistory,
2766                      gmx_bool                       reproducibilityRequested,
2767                      const gmx::MdModulesNotifier&  mdModulesNotifier,
2768                      gmx::ReadCheckpointDataHolder* modularSimulatorCheckpointData,
2769                      bool                           useModularSimulator)
2770 {
2771     CheckpointHeaderContents headerContents;
2772     if (SIMMASTER(cr))
2773     {
2774         /* Read the state from the checkpoint file */
2775         read_checkpoint(fn,
2776                         logfio,
2777                         cr,
2778                         dd_nc,
2779                         ir->eI,
2780                         &(ir->fepvals->init_fep_state),
2781                         &headerContents,
2782                         state,
2783                         observablesHistory,
2784                         reproducibilityRequested,
2785                         mdModulesNotifier,
2786                         modularSimulatorCheckpointData,
2787                         useModularSimulator);
2788     }
2789     if (PAR(cr))
2790     {
2791         gmx_bcast(sizeof(headerContents.step), &headerContents.step, cr->mpiDefaultCommunicator);
2792         gmx::MdModulesCheckpointReadingBroadcast broadcastCheckPointData = {
2793             cr->mpiDefaultCommunicator, PAR(cr), headerContents.file_version
2794         };
2795         mdModulesNotifier.checkpointingNotifications_.notify(broadcastCheckPointData);
2796     }
2797     ir->bContinuation = TRUE;
2798     if (ir->nsteps >= 0)
2799     {
2800         // TODO Should the following condition be <=? Currently if you
2801         // pass a checkpoint written by an normal completion to a restart,
2802         // mdrun will read all input, does some work but no steps, and
2803         // write successful output. But perhaps that is not desirable.
2804         // Note that we do not intend to support the use of mdrun
2805         // -nsteps to circumvent this condition.
2806         if (ir->nsteps + ir->init_step < headerContents.step)
2807         {
2808             char        buf[STEPSTRSIZE];
2809             std::string message =
2810                     gmx::formatString("The input requested %s steps, ", gmx_step_str(ir->nsteps, buf));
2811             if (ir->init_step > 0)
2812             {
2813                 message += gmx::formatString("starting from step %s, ", gmx_step_str(ir->init_step, buf));
2814             }
2815             message += gmx::formatString(
2816                     "however the checkpoint "
2817                     "file has already reached step %s. The simulation will not "
2818                     "proceed, because either your simulation is already complete, "
2819                     "or your combination of input files don't match.",
2820                     gmx_step_str(headerContents.step, buf));
2821             gmx_fatal(FARGS, "%s", message.c_str());
2822         }
2823         ir->nsteps += ir->init_step - headerContents.step;
2824     }
2825     ir->init_step       = headerContents.step;
2826     ir->simulation_part = headerContents.simulation_part + 1;
2827 }
2828
2829 void read_checkpoint_part_and_step(const char* filename, int* simulation_part, int64_t* step)
2830 {
2831     t_fileio* fp;
2832
2833     if (filename == nullptr || !gmx_fexist(filename) || ((fp = gmx_fio_open(filename, "r")) == nullptr))
2834     {
2835         *simulation_part = 0;
2836         *step            = 0;
2837         return;
2838     }
2839
2840     CheckpointHeaderContents headerContents;
2841     do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, &headerContents);
2842     gmx_fio_close(fp);
2843     *simulation_part = headerContents.simulation_part;
2844     *step            = headerContents.step;
2845 }
2846
2847 static CheckpointHeaderContents read_checkpoint_data(t_fileio*                         fp,
2848                                                      t_state*                          state,
2849                                                      std::vector<gmx_file_position_t>* outputfiles,
2850                                                      gmx::ReadCheckpointDataHolder* modularSimulatorCheckpointData)
2851 {
2852     CheckpointHeaderContents headerContents;
2853     do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, &headerContents);
2854     state->natoms        = headerContents.natoms;
2855     state->ngtc          = headerContents.ngtc;
2856     state->nnhpres       = headerContents.nnhpres;
2857     state->nhchainlength = headerContents.nhchainlength;
2858     state->flags         = headerContents.flags_state;
2859     int ret              = do_cpt_state(gmx_fio_getxdr(fp), state->flags, state, nullptr);
2860     if (ret)
2861     {
2862         cp_error();
2863     }
2864     ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents.flags_eks, &state->ekinstate, nullptr);
2865     if (ret)
2866     {
2867         cp_error();
2868     }
2869
2870     energyhistory_t enerhist;
2871     ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE, headerContents.flags_enh, &enerhist, nullptr);
2872     if (ret)
2873     {
2874         cp_error();
2875     }
2876     PullHistory pullHist = {};
2877     ret                  = doCptPullHist(
2878             gmx_fio_getxdr(fp), TRUE, headerContents.flagsPullHistory, &pullHist, StatePart::pullHistory, nullptr);
2879     if (ret)
2880     {
2881         cp_error();
2882     }
2883
2884     ret = do_cpt_df_hist(
2885             gmx_fio_getxdr(fp), headerContents.flags_dfh, headerContents.nlambda, &state->dfhist, nullptr);
2886     if (ret)
2887     {
2888         cp_error();
2889     }
2890
2891     edsamhistory_t edsamhist = {};
2892     ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, headerContents.nED, &edsamhist, nullptr);
2893     if (ret)
2894     {
2895         cp_error();
2896     }
2897
2898     ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE, headerContents.flags_awhh, state->awhHistory.get(), nullptr);
2899     if (ret)
2900     {
2901         cp_error();
2902     }
2903
2904     swaphistory_t swaphist = {};
2905     ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents.eSwapCoords, &swaphist, nullptr);
2906     if (ret)
2907     {
2908         cp_error();
2909     }
2910
2911     ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, outputfiles, nullptr, headerContents.file_version);
2912
2913     if (ret)
2914     {
2915         cp_error();
2916     }
2917     gmx::MdModulesNotifier mdModuleNotifier;
2918     do_cpt_mdmodules(headerContents.file_version, fp, mdModuleNotifier, nullptr);
2919     if (headerContents.file_version >= cptv_ModularSimulator)
2920     {
2921         // Store modular checkpoint data into modularSimulatorCheckpointData
2922         gmx::FileIOXdrSerializer serializer(fp);
2923         modularSimulatorCheckpointData->deserialize(&serializer);
2924     }
2925     ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
2926     if (ret)
2927     {
2928         cp_error();
2929     }
2930     return headerContents;
2931 }
2932
2933 void read_checkpoint_trxframe(t_fileio* fp, t_trxframe* fr)
2934 {
2935     t_state                          state;
2936     std::vector<gmx_file_position_t> outputfiles;
2937     gmx::ReadCheckpointDataHolder    modularSimulatorCheckpointData;
2938     CheckpointHeaderContents         headerContents =
2939             read_checkpoint_data(fp, &state, &outputfiles, &modularSimulatorCheckpointData);
2940     if (headerContents.isModularSimulatorCheckpoint)
2941     {
2942         gmx::ModularSimulator::readCheckpointToTrxFrame(fr, &modularSimulatorCheckpointData, headerContents);
2943         return;
2944     }
2945
2946     fr->natoms    = state.natoms;
2947     fr->bStep     = TRUE;
2948     fr->step      = int64_to_int(headerContents.step, "conversion of checkpoint to trajectory");
2949     fr->bTime     = TRUE;
2950     fr->time      = headerContents.t;
2951     fr->bLambda   = TRUE;
2952     fr->lambda    = state.lambda[efptFEP];
2953     fr->fep_state = state.fep_state;
2954     fr->bAtoms    = FALSE;
2955     fr->bX        = ((state.flags & (1 << estX)) != 0);
2956     if (fr->bX)
2957     {
2958         fr->x = makeRvecArray(state.x, state.natoms);
2959     }
2960     fr->bV = ((state.flags & (1 << estV)) != 0);
2961     if (fr->bV)
2962     {
2963         fr->v = makeRvecArray(state.v, state.natoms);
2964     }
2965     fr->bF   = FALSE;
2966     fr->bBox = ((state.flags & (1 << estBOX)) != 0);
2967     if (fr->bBox)
2968     {
2969         copy_mat(state.box, fr->box);
2970     }
2971 }
2972
2973 void list_checkpoint(const char* fn, FILE* out)
2974 {
2975     t_fileio* fp;
2976     int       ret;
2977
2978     t_state state;
2979
2980     fp = gmx_fio_open(fn, "r");
2981     CheckpointHeaderContents headerContents;
2982     do_cpt_header(gmx_fio_getxdr(fp), TRUE, out, &headerContents);
2983     state.natoms        = headerContents.natoms;
2984     state.ngtc          = headerContents.ngtc;
2985     state.nnhpres       = headerContents.nnhpres;
2986     state.nhchainlength = headerContents.nhchainlength;
2987     state.flags         = headerContents.flags_state;
2988     ret                 = do_cpt_state(gmx_fio_getxdr(fp), state.flags, &state, out);
2989     if (ret)
2990     {
2991         cp_error();
2992     }
2993     ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents.flags_eks, &state.ekinstate, out);
2994     if (ret)
2995     {
2996         cp_error();
2997     }
2998
2999     energyhistory_t enerhist;
3000     ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE, headerContents.flags_enh, &enerhist, out);
3001
3002     if (ret == 0)
3003     {
3004         PullHistory pullHist = {};
3005         ret                  = doCptPullHist(
3006                 gmx_fio_getxdr(fp), TRUE, headerContents.flagsPullHistory, &pullHist, StatePart::pullHistory, out);
3007     }
3008
3009     if (ret == 0)
3010     {
3011         ret = do_cpt_df_hist(
3012                 gmx_fio_getxdr(fp), headerContents.flags_dfh, headerContents.nlambda, &state.dfhist, out);
3013     }
3014
3015     if (ret == 0)
3016     {
3017         edsamhistory_t edsamhist = {};
3018         ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, headerContents.nED, &edsamhist, out);
3019     }
3020
3021     if (ret == 0)
3022     {
3023         ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE, headerContents.flags_awhh, state.awhHistory.get(), out);
3024     }
3025
3026     if (ret == 0)
3027     {
3028         swaphistory_t swaphist = {};
3029         ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents.eSwapCoords, &swaphist, out);
3030     }
3031
3032     if (ret == 0)
3033     {
3034         std::vector<gmx_file_position_t> outputfiles;
3035         ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, out, headerContents.file_version);
3036     }
3037     gmx::MdModulesNotifier mdModuleNotifier;
3038     do_cpt_mdmodules(headerContents.file_version, fp, mdModuleNotifier, out);
3039     if (headerContents.file_version >= cptv_ModularSimulator)
3040     {
3041         gmx::FileIOXdrSerializer      serializer(fp);
3042         gmx::ReadCheckpointDataHolder modularSimulatorCheckpointData;
3043         modularSimulatorCheckpointData.deserialize(&serializer);
3044         modularSimulatorCheckpointData.dump(out);
3045     }
3046
3047     if (ret == 0)
3048     {
3049         ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
3050     }
3051
3052     if (ret)
3053     {
3054         cp_warning(out);
3055     }
3056     if (gmx_fio_close(fp) != 0)
3057     {
3058         gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
3059     }
3060 }
3061
3062 /* This routine cannot print tons of data, since it is called before the log file is opened. */
3063 CheckpointHeaderContents read_checkpoint_simulation_part_and_filenames(t_fileio* fp,
3064                                                                        std::vector<gmx_file_position_t>* outputfiles)
3065 {
3066     t_state                       state;
3067     gmx::ReadCheckpointDataHolder modularSimulatorCheckpointData;
3068     CheckpointHeaderContents      headerContents =
3069             read_checkpoint_data(fp, &state, outputfiles, &modularSimulatorCheckpointData);
3070     if (gmx_fio_close(fp) != 0)
3071     {
3072         gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
3073     }
3074     return headerContents;
3075 }