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