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