src/gromacs/fileio/checkpoint.cpp

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