src/gromacs/fileio/checkpoint.cpp

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