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