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