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     std::unique_ptr<df_history_t> localDFHistory = nullptr;
1908     if (*dfhistPtr == nullptr)
1909     {
1910         localDFHistory        = std::make_unique<df_history_t>();
1911         *dfhistPtr            = localDFHistory.get();
1912         (*dfhistPtr)->nlambda = nlambda;
1913         init_df_history(*dfhistPtr, nlambda);
1914     }
1915     df_history_t* dfhist = *dfhistPtr;
1916
1917     using StateFlags = gmx::EnumerationArray<StateFepEntry, bool>;
1918     for (auto i = StateFlags::keys().begin(); (i != StateFlags::keys().end() && ret == 0); i++)
1919     {
1920         if (fflags & enumValueToBitMask(*i))
1921         {
1922             switch (*i)
1923             {
1924                 case StateFepEntry::IsEquilibrated:
1925                     ret = do_cpte_bool(xd, *i, fflags, &dfhist->bEquil, list);
1926                     break;
1927                 case StateFepEntry::NumAtLambda:
1928                     ret = do_cpte_ints(xd, *i, fflags, nlambda, &dfhist->n_at_lam, list);
1929                     break;
1930                 case StateFepEntry::WangLandauHistogram:
1931                     ret = do_cpte_reals(xd, *i, fflags, nlambda, &dfhist->wl_histo, list);
1932                     break;
1933                 case StateFepEntry::WangLandauDelta:
1934                     ret = do_cpte_real(xd, *i, fflags, &dfhist->wl_delta, list);
1935                     break;
1936                 case StateFepEntry::SumWeights:
1937                     ret = do_cpte_reals(xd, *i, fflags, nlambda, &dfhist->sum_weights, list);
1938                     break;
1939                 case StateFepEntry::SumDG:
1940                     ret = do_cpte_reals(xd, *i, fflags, nlambda, &dfhist->sum_dg, list);
1941                     break;
1942                 case StateFepEntry::SumMinVar:
1943                     ret = do_cpte_reals(xd, *i, fflags, nlambda, &dfhist->sum_minvar, list);
1944                     break;
1945                 case StateFepEntry::SumVar:
1946                     ret = do_cpte_reals(xd, *i, fflags, nlambda, &dfhist->sum_variance, list);
1947                     break;
1948                 case StateFepEntry::Accump:
1949                     ret = do_cpte_nmatrix(xd, *i, fflags, nlambda, dfhist->accum_p, list);
1950                     break;
1951                 case StateFepEntry::Accumm:
1952                     ret = do_cpte_nmatrix(xd, *i, fflags, nlambda, dfhist->accum_m, list);
1953                     break;
1954                 case StateFepEntry::Accump2:
1955                     ret = do_cpte_nmatrix(xd, *i, fflags, nlambda, dfhist->accum_p2, list);
1956                     break;
1957                 case StateFepEntry::Accumm2:
1958                     ret = do_cpte_nmatrix(xd, *i, fflags, nlambda, dfhist->accum_m2, list);
1959                     break;
1960                 case StateFepEntry::Tij:
1961                     ret = do_cpte_nmatrix(xd, *i, fflags, nlambda, dfhist->Tij, list);
1962                     break;
1963                 case StateFepEntry::TijEmp:
1964                     ret = do_cpte_nmatrix(xd, *i, fflags, nlambda, dfhist->Tij_empirical, list);
1965                     break;
1966
1967                 default:
1968                     gmx_fatal(FARGS,
1969                               "Unknown df history entry %d\n"
1970                               "You are probably reading a new checkpoint file with old code",
1971                               enumValueToBitMask(*i));
1972             }
1973         }
1974     }
1975
1976     return ret;
1977 }
1978
1979
1980 /* This function stores the last whole configuration of the reference and
1981  * average structure in the .cpt file
1982  */
1983 static int do_cpt_EDstate(XDR* xd, gmx_bool bRead, int nED, edsamhistory_t* EDstate, FILE* list)
1984 {
1985     if (nED == 0)
1986     {
1987         return 0;
1988     }
1989
1990     EDstate->bFromCpt = bRead;
1991     EDstate->nED      = nED;
1992
1993     /* When reading, init_edsam has not been called yet,
1994      * so we have to allocate memory first. */
1995     if (bRead)
1996     {
1997         snew(EDstate->nref, EDstate->nED);
1998         snew(EDstate->old_sref, EDstate->nED);
1999         snew(EDstate->nav, EDstate->nED);
2000         snew(EDstate->old_sav, EDstate->nED);
2001     }
2002
2003     /* Read/write the last whole conformation of SREF and SAV for each ED dataset (usually only one) */
2004     for (int i = 0; i < EDstate->nED; i++)
2005     {
2006         char buf[STRLEN];
2007
2008         /* Reference structure SREF */
2009         sprintf(buf, "ED%d # of atoms in reference structure", i + 1);
2010         do_cpt_int_err(xd, buf, &EDstate->nref[i], list);
2011         sprintf(buf, "ED%d x_ref", i + 1);
2012         if (bRead)
2013         {
2014             snew(EDstate->old_sref[i], EDstate->nref[i]);
2015             do_cpt_n_rvecs_err(xd, buf, EDstate->nref[i], EDstate->old_sref[i], list);
2016         }
2017         else
2018         {
2019             do_cpt_n_rvecs_err(xd, buf, EDstate->nref[i], EDstate->old_sref_p[i], list);
2020         }
2021
2022         /* Average structure SAV */
2023         sprintf(buf, "ED%d # of atoms in average structure", i + 1);
2024         do_cpt_int_err(xd, buf, &EDstate->nav[i], list);
2025         sprintf(buf, "ED%d x_av", i + 1);
2026         if (bRead)
2027         {
2028             snew(EDstate->old_sav[i], EDstate->nav[i]);
2029             do_cpt_n_rvecs_err(xd, buf, EDstate->nav[i], EDstate->old_sav[i], list);
2030         }
2031         else
2032         {
2033             do_cpt_n_rvecs_err(xd, buf, EDstate->nav[i], EDstate->old_sav_p[i], list);
2034         }
2035     }
2036
2037     return 0;
2038 }
2039
2040 static int do_cpt_correlation_grid(XDR*                         xd,
2041                                    gmx_bool                     bRead,
2042                                    gmx_unused int               fflags,
2043                                    gmx::CorrelationGridHistory* corrGrid,
2044                                    FILE*                        list,
2045                                    StateAwhEntry                eawhh)
2046 {
2047     int ret = 0;
2048
2049     do_cpt_int_err(xd, enumValueToString(eawhh), &(corrGrid->numCorrelationTensors), list);
2050     do_cpt_int_err(xd, enumValueToString(eawhh), &(corrGrid->tensorSize), list);
2051     do_cpt_int_err(xd, enumValueToString(eawhh), &(corrGrid->blockDataListSize), list);
2052
2053     if (bRead)
2054     {
2055         initCorrelationGridHistory(
2056                 corrGrid, corrGrid->numCorrelationTensors, corrGrid->tensorSize, corrGrid->blockDataListSize);
2057     }
2058
2059     for (gmx::CorrelationBlockDataHistory& blockData : corrGrid->blockDataBuffer)
2060     {
2061         do_cpt_double_err(xd, enumValueToString(eawhh), &(blockData.blockSumWeight), list);
2062         do_cpt_double_err(xd, enumValueToString(eawhh), &(blockData.blockSumSquareWeight), list);
2063         do_cpt_double_err(xd, enumValueToString(eawhh), &(blockData.blockSumWeightX), list);
2064         do_cpt_double_err(xd, enumValueToString(eawhh), &(blockData.blockSumWeightY), list);
2065         do_cpt_double_err(xd, enumValueToString(eawhh), &(blockData.sumOverBlocksSquareBlockWeight), list);
2066         do_cpt_double_err(xd, enumValueToString(eawhh), &(blockData.sumOverBlocksBlockSquareWeight), list);
2067         do_cpt_double_err(
2068                 xd, enumValueToString(eawhh), &(blockData.sumOverBlocksBlockWeightBlockWeightX), list);
2069         do_cpt_double_err(
2070                 xd, enumValueToString(eawhh), &(blockData.sumOverBlocksBlockWeightBlockWeightY), list);
2071         do_cpt_double_err(xd, enumValueToString(eawhh), &(blockData.blockLength), list);
2072         do_cpt_int_err(xd, enumValueToString(eawhh), &(blockData.previousBlockIndex), list);
2073         do_cpt_double_err(xd, enumValueToString(eawhh), &(blockData.correlationIntegral), list);
2074     }
2075
2076     return ret;
2077 }
2078
2079 static int do_cpt_awh_bias(XDR* xd, gmx_bool bRead, int fflags, gmx::AwhBiasHistory* biasHistory, FILE* list)
2080 {
2081     int ret = 0;
2082
2083     gmx::AwhBiasStateHistory* state = &biasHistory->state;
2084     using StateFlags                = gmx::EnumerationArray<StateAwhEntry, bool>;
2085     for (auto i = StateFlags::keys().begin(); (i != StateFlags::keys().end() && ret == 0); i++)
2086     {
2087         if (fflags & enumValueToBitMask(*i))
2088         {
2089             switch (*i)
2090             {
2091                 case StateAwhEntry::InInitial:
2092                     do_cpt_bool_err(xd, enumValueToString(*i), &state->in_initial, list);
2093                     break;
2094                 case StateAwhEntry::EquilibrateHistogram:
2095                     do_cpt_bool_err(xd, enumValueToString(*i), &state->equilibrateHistogram, list);
2096                     break;
2097                 case StateAwhEntry::HistogramSize:
2098                     do_cpt_double_err(xd, enumValueToString(*i), &state->histSize, list);
2099                     break;
2100                 case StateAwhEntry::NumPoints:
2101                 {
2102                     int numPoints;
2103                     if (!bRead)
2104                     {
2105                         numPoints = biasHistory->pointState.size();
2106                     }
2107                     do_cpt_int_err(xd, enumValueToString(*i), &numPoints, list);
2108                     if (bRead)
2109                     {
2110                         biasHistory->pointState.resize(numPoints);
2111                     }
2112                 }
2113                 break;
2114                 case StateAwhEntry::CoordPoint:
2115                     for (auto& psh : biasHistory->pointState)
2116                     {
2117                         do_cpt_double_err(xd, enumValueToString(*i), &psh.target, list);
2118                         do_cpt_double_err(xd, enumValueToString(*i), &psh.free_energy, list);
2119                         do_cpt_double_err(xd, enumValueToString(*i), &psh.bias, list);
2120                         do_cpt_double_err(xd, enumValueToString(*i), &psh.weightsum_iteration, list);
2121                         do_cpt_double_err(xd, enumValueToString(*i), &psh.weightsum_covering, list);
2122                         do_cpt_double_err(xd, enumValueToString(*i), &psh.weightsum_tot, list);
2123                         do_cpt_double_err(xd, enumValueToString(*i), &psh.weightsum_ref, list);
2124                         do_cpt_step_err(xd, enumValueToString(*i), &psh.last_update_index, list);
2125                         do_cpt_double_err(xd, enumValueToString(*i), &psh.log_pmfsum, list);
2126                         do_cpt_double_err(xd, enumValueToString(*i), &psh.visits_iteration, list);
2127                         do_cpt_double_err(xd, enumValueToString(*i), &psh.visits_tot, list);
2128                     }
2129                     break;
2130                 case StateAwhEntry::UmbrellaGridPoint:
2131                     do_cpt_int_err(xd, enumValueToString(*i), &(state->umbrellaGridpoint), list);
2132                     break;
2133                 case StateAwhEntry::UpdateList:
2134                     do_cpt_int_err(xd, enumValueToString(*i), &(state->origin_index_updatelist), list);
2135                     do_cpt_int_err(xd, enumValueToString(*i), &(state->end_index_updatelist), list);
2136                     break;
2137                 case StateAwhEntry::LogScaledSampleWeight:
2138                     do_cpt_double_err(xd, enumValueToString(*i), &(state->logScaledSampleWeight), list);
2139                     do_cpt_double_err(xd, enumValueToString(*i), &(state->maxLogScaledSampleWeight), list);
2140                     break;
2141                 case StateAwhEntry::NumUpdates:
2142                     do_cpt_step_err(xd, enumValueToString(*i), &(state->numUpdates), list);
2143                     break;
2144                 case StateAwhEntry::ForceCorrelationGrid:
2145                     ret = do_cpt_correlation_grid(
2146                             xd, bRead, fflags, &biasHistory->forceCorrelationGrid, list, *i);
2147                     break;
2148                 default: gmx_fatal(FARGS, "Unknown awh history entry %d\n", enumValueToBitMask(*i));
2149             }
2150         }
2151     }
2152
2153     return ret;
2154 }
2155
2156 static int do_cpt_awh(XDR* xd, gmx_bool bRead, int fflags, gmx::AwhHistory* awhHistory, FILE* list)
2157 {
2158     int ret = 0;
2159
2160     if (fflags != 0)
2161     {
2162         std::shared_ptr<gmx::AwhHistory> awhHistoryLocal;
2163
2164         if (awhHistory == nullptr)
2165         {
2166             GMX_RELEASE_ASSERT(bRead,
2167                                "do_cpt_awh should not be called for writing without an AwhHistory");
2168
2169             awhHistoryLocal = std::make_shared<gmx::AwhHistory>();
2170             awhHistory      = awhHistoryLocal.get();
2171         }
2172
2173         /* To be able to read and write the AWH data several parameters determining the layout of the AWH structs need to be known
2174            (nbias, npoints, etc.). The best thing (?) would be to have these passed to this function. When writing to a checkpoint
2175            these parameters are available in awh_history (after calling init_awh_history). When reading from a checkpoint though, there
2176            is no initialized awh_history (it is initialized and set in this function). The AWH parameters have not always been read
2177            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
2178            when writing a checkpoint, also storing parameters needed for future reading of the checkpoint.
2179
2180            Another issue is that some variables that AWH checkpoints don't have a registered enum and string (e.g. nbias below).
2181            One difficulty is the multilevel structure of the data which would need to be represented somehow. */
2182
2183         int numBias;
2184         if (!bRead)
2185         {
2186             numBias = awhHistory->bias.size();
2187         }
2188         do_cpt_int_err(xd, "awh_nbias", &numBias, list);
2189
2190         if (bRead)
2191         {
2192             awhHistory->bias.resize(numBias);
2193         }
2194         for (auto& bias : awhHistory->bias)
2195         {
2196             ret = do_cpt_awh_bias(xd, bRead, fflags, &bias, list);
2197             if (ret)
2198             {
2199                 return ret;
2200             }
2201         }
2202         do_cpt_double_err(xd, "awh_potential_offset", &awhHistory->potentialOffset, list);
2203     }
2204     return ret;
2205 }
2206
2207 static void do_cpt_mdmodules(CheckPointVersion              fileVersion,
2208                              t_fileio*                      checkpointFileHandle,
2209                              const gmx::MDModulesNotifiers& mdModulesNotifiers,
2210                              FILE*                          outputFile)
2211 {
2212     if (fileVersion >= CheckPointVersion::MDModules)
2213     {
2214         gmx::FileIOXdrSerializer serializer(checkpointFileHandle);
2215         gmx::KeyValueTreeObject  mdModuleCheckpointParameterTree =
2216                 gmx::deserializeKeyValueTree(&serializer);
2217         if (outputFile)
2218         {
2219             gmx::TextWriter textWriter(outputFile);
2220             gmx::dumpKeyValueTree(&textWriter, mdModuleCheckpointParameterTree);
2221         }
2222         gmx::MDModulesCheckpointReadingDataOnMaster mdModuleCheckpointReadingDataOnMaster = {
2223             mdModuleCheckpointParameterTree, fileVersion
2224         };
2225         mdModulesNotifiers.checkpointingNotifier_.notify(mdModuleCheckpointReadingDataOnMaster);
2226     }
2227 }
2228
2229 static int do_cpt_files(XDR*                              xd,
2230                         gmx_bool                          bRead,
2231                         std::vector<gmx_file_position_t>* outputfiles,
2232                         FILE*                             list,
2233                         CheckPointVersion                 file_version)
2234 {
2235     gmx_off_t                   offset;
2236     gmx_off_t                   mask = 0xFFFFFFFFL;
2237     int                         offset_high, offset_low;
2238     std::array<char, CPTSTRLEN> buf;
2239     GMX_RELEASE_ASSERT(outputfiles, "Must have valid outputfiles");
2240
2241     // Ensure that reading pre-allocates outputfiles, while writing
2242     // writes what is already there.
2243     int nfiles = outputfiles->size();
2244     if (do_cpt_int(xd, "number of output files", &nfiles, list) != 0)
2245     {
2246         return -1;
2247     }
2248     if (bRead)
2249     {
2250         outputfiles->resize(nfiles);
2251     }
2252
2253     for (auto& outputfile : *outputfiles)
2254     {
2255         /* 64-bit XDR numbers are not portable, so it is stored as separate high/low fractions */
2256         if (bRead)
2257         {
2258             do_cpt_string_err(xd, "output filename", buf, list);
2259             std::copy(std::begin(buf), std::end(buf), std::begin(outputfile.filename));
2260
2261             if (do_cpt_int(xd, "file_offset_high", &offset_high, list) != 0)
2262             {
2263                 return -1;
2264             }
2265             if (do_cpt_int(xd, "file_offset_low", &offset_low, list) != 0)
2266             {
2267                 return -1;
2268             }
2269             outputfile.offset = (static_cast<gmx_off_t>(offset_high) << 32)
2270                                 | (static_cast<gmx_off_t>(offset_low) & mask);
2271         }
2272         else
2273         {
2274             do_cpt_string_err(xd, "output filename", outputfile.filename, list);
2275             /* writing */
2276             offset = outputfile.offset;
2277             if (offset == -1)
2278             {
2279                 offset_low  = -1;
2280                 offset_high = -1;
2281             }
2282             else
2283             {
2284                 offset_low  = static_cast<int>(offset & mask);
2285                 offset_high = static_cast<int>((offset >> 32) & mask);
2286             }
2287             if (do_cpt_int(xd, "file_offset_high", &offset_high, list) != 0)
2288             {
2289                 return -1;
2290             }
2291             if (do_cpt_int(xd, "file_offset_low", &offset_low, list) != 0)
2292             {
2293                 return -1;
2294             }
2295         }
2296         if (file_version >= CheckPointVersion::FileChecksumAndSize)
2297         {
2298             if (do_cpt_int(xd, "file_checksum_size", &outputfile.checksumSize, list) != 0)
2299             {
2300                 return -1;
2301             }
2302             if (do_cpt_u_chars(xd, "file_checksum", outputfile.checksum.size(), outputfile.checksum.data(), list)
2303                 != 0)
2304             {
2305                 return -1;
2306             }
2307         }
2308         else
2309         {
2310             outputfile.checksumSize = -1;
2311         }
2312     }
2313     return 0;
2314 }
2315
2316 void write_checkpoint_data(t_fileio*                         fp,
2317                            CheckpointHeaderContents          headerContents,
2318                            gmx_bool                          bExpanded,
2319                            LambdaWeightCalculation           elamstats,
2320                            t_state*                          state,
2321                            ObservablesHistory*               observablesHistory,
2322                            const gmx::MDModulesNotifiers&    mdModulesNotifiers,
2323                            std::vector<gmx_file_position_t>* outputfiles,
2324                            gmx::WriteCheckpointDataHolder*   modularSimulatorCheckpointData)
2325 {
2326     headerContents.flags_eks = 0;
2327     if (state->ekinstate.bUpToDate)
2328     {
2329         headerContents.flags_eks = (enumValueToBitMask(StateKineticEntry::EkinNumber)
2330                                     | enumValueToBitMask(StateKineticEntry::EkinHalfStep)
2331                                     | enumValueToBitMask(StateKineticEntry::EkinFullStep)
2332                                     | enumValueToBitMask(StateKineticEntry::EkinHalfStepOld)
2333                                     | enumValueToBitMask(StateKineticEntry::EkinNoseHooverScaleFullStep)
2334                                     | enumValueToBitMask(StateKineticEntry::EkinNoseHooverScaleHalfStep)
2335                                     | enumValueToBitMask(StateKineticEntry::VelocityScale)
2336                                     | enumValueToBitMask(StateKineticEntry::DEkinDLambda)
2337                                     | enumValueToBitMask(StateKineticEntry::Mvcos));
2338     }
2339     headerContents.isModularSimulatorCheckpoint = !modularSimulatorCheckpointData->empty();
2340
2341     energyhistory_t* enerhist = observablesHistory->energyHistory.get();
2342     headerContents.flags_enh  = 0;
2343     if (enerhist != nullptr && (enerhist->nsum > 0 || enerhist->nsum_sim > 0))
2344     {
2345         headerContents.flags_enh |= enumValueToBitMask(StateEnergyEntry::N)
2346                                     | enumValueToBitMask(StateEnergyEntry::NumSteps)
2347                                     | enumValueToBitMask(StateEnergyEntry::NumStepsSim);
2348         if (enerhist->nsum > 0)
2349         {
2350             headerContents.flags_enh |= (enumValueToBitMask(StateEnergyEntry::Aver)
2351                                          | enumValueToBitMask(StateEnergyEntry::Sum)
2352                                          | enumValueToBitMask(StateEnergyEntry::NumSum));
2353         }
2354         if (enerhist->nsum_sim > 0)
2355         {
2356             headerContents.flags_enh |= (enumValueToBitMask(StateEnergyEntry::SumSim)
2357                                          | enumValueToBitMask(StateEnergyEntry::NumSumSim));
2358         }
2359         if (enerhist->deltaHForeignLambdas != nullptr)
2360         {
2361             headerContents.flags_enh |= (enumValueToBitMask(StateEnergyEntry::DeltaHNN)
2362                                          | enumValueToBitMask(StateEnergyEntry::DeltaHList)
2363                                          | enumValueToBitMask(StateEnergyEntry::DeltaHStartTime)
2364                                          | enumValueToBitMask(StateEnergyEntry::DeltaHStartLambda));
2365         }
2366     }
2367
2368     PullHistory* pullHist           = observablesHistory->pullHistory.get();
2369     headerContents.flagsPullHistory = 0;
2370     if (pullHist != nullptr && (pullHist->numValuesInXSum > 0 || pullHist->numValuesInFSum > 0))
2371     {
2372         headerContents.flagsPullHistory |= enumValueToBitMask(StatePullEntry::NumCoordinates);
2373         headerContents.flagsPullHistory |= (enumValueToBitMask(StatePullEntry::NumGroups)
2374                                             | enumValueToBitMask(StatePullEntry::NumValuesInXSum)
2375                                             | enumValueToBitMask(StatePullEntry::NumValuesInFSum));
2376     }
2377
2378     headerContents.flags_dfh = 0;
2379     if (bExpanded)
2380     {
2381         headerContents.flags_dfh =
2382                 (enumValueToBitMask(StateFepEntry::IsEquilibrated)
2383                  | enumValueToBitMask(StateFepEntry::NumAtLambda)
2384                  | enumValueToBitMask(StateFepEntry::SumWeights) | enumValueToBitMask(StateFepEntry::SumDG)
2385                  | enumValueToBitMask(StateFepEntry::Tij) | enumValueToBitMask(StateFepEntry::TijEmp));
2386         if (EWL(elamstats))
2387         {
2388             headerContents.flags_dfh |= (enumValueToBitMask(StateFepEntry::WangLandauDelta)
2389                                          | enumValueToBitMask(StateFepEntry::WangLandauHistogram));
2390         }
2391         if ((elamstats == LambdaWeightCalculation::Minvar) || (elamstats == LambdaWeightCalculation::Barker)
2392             || (elamstats == LambdaWeightCalculation::Metropolis))
2393         {
2394             headerContents.flags_dfh |= (enumValueToBitMask(StateFepEntry::Accump)
2395                                          | enumValueToBitMask(StateFepEntry::Accumm)
2396                                          | enumValueToBitMask(StateFepEntry::Accump2)
2397                                          | enumValueToBitMask(StateFepEntry::Accumm2)
2398                                          | enumValueToBitMask(StateFepEntry::SumMinVar)
2399                                          | enumValueToBitMask(StateFepEntry::SumVar));
2400         }
2401     }
2402
2403     headerContents.flags_awhh = 0;
2404     if (state->awhHistory != nullptr && !state->awhHistory->bias.empty())
2405     {
2406         headerContents.flags_awhh |= (enumValueToBitMask(StateAwhEntry::InInitial)
2407                                       | enumValueToBitMask(StateAwhEntry::EquilibrateHistogram)
2408                                       | enumValueToBitMask(StateAwhEntry::HistogramSize)
2409                                       | enumValueToBitMask(StateAwhEntry::NumPoints)
2410                                       | enumValueToBitMask(StateAwhEntry::CoordPoint)
2411                                       | enumValueToBitMask(StateAwhEntry::UmbrellaGridPoint)
2412                                       | enumValueToBitMask(StateAwhEntry::UpdateList)
2413                                       | enumValueToBitMask(StateAwhEntry::LogScaledSampleWeight)
2414                                       | enumValueToBitMask(StateAwhEntry::NumUpdates)
2415                                       | enumValueToBitMask(StateAwhEntry::ForceCorrelationGrid));
2416     }
2417
2418     do_cpt_header(gmx_fio_getxdr(fp), FALSE, nullptr, &headerContents);
2419
2420     if ((do_cpt_state(gmx_fio_getxdr(fp), state->flags, state, nullptr) < 0)
2421         || (do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents.flags_eks, &state->ekinstate, nullptr) < 0)
2422         || (do_cpt_enerhist(gmx_fio_getxdr(fp), FALSE, headerContents.flags_enh, enerhist, nullptr) < 0)
2423         || (doCptPullHist(gmx_fio_getxdr(fp), FALSE, headerContents.flagsPullHistory, pullHist, nullptr) < 0)
2424         || (do_cpt_df_hist(gmx_fio_getxdr(fp), headerContents.flags_dfh, headerContents.nlambda, &state->dfhist, nullptr)
2425             < 0)
2426         || (do_cpt_EDstate(
2427                     gmx_fio_getxdr(fp), FALSE, headerContents.nED, observablesHistory->edsamHistory.get(), nullptr)
2428             < 0)
2429         || (do_cpt_awh(gmx_fio_getxdr(fp), FALSE, headerContents.flags_awhh, state->awhHistory.get(), nullptr) < 0)
2430         || (do_cpt_swapstate(gmx_fio_getxdr(fp),
2431                              FALSE,
2432                              headerContents.eSwapCoords,
2433                              observablesHistory->swapHistory.get(),
2434                              nullptr)
2435             < 0)
2436         || (do_cpt_files(gmx_fio_getxdr(fp), FALSE, outputfiles, nullptr, headerContents.file_version) < 0))
2437     {
2438         gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
2439     }
2440
2441     // Checkpointing MDModules
2442     {
2443         gmx::KeyValueTreeBuilder          builder;
2444         gmx::MDModulesWriteCheckpointData mdModulesWriteCheckpoint = { builder.rootObject(),
2445                                                                        headerContents.file_version };
2446         mdModulesNotifiers.checkpointingNotifier_.notify(mdModulesWriteCheckpoint);
2447         auto                     tree = builder.build();
2448         gmx::FileIOXdrSerializer serializer(fp);
2449         gmx::serializeKeyValueTree(tree, &serializer);
2450     }
2451
2452     // Checkpointing modular simulator
2453     {
2454         gmx::FileIOXdrSerializer serializer(fp);
2455         modularSimulatorCheckpointData->serialize(&serializer);
2456     }
2457
2458     do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
2459 #if GMX_FAHCORE
2460     /* Always FAH checkpoint immediately after a Gromacs checkpoint.
2461      *
2462      * Note that it is critical that we save a FAH checkpoint directly
2463      * after writing a Gromacs checkpoint.  If the program dies, either
2464      * by the machine powering off suddenly or the process being,
2465      * killed, FAH can recover files that have only appended data by
2466      * truncating them to the last recorded length.  The Gromacs
2467      * checkpoint does not just append data, it is fully rewritten each
2468      * time so a crash between moving the new Gromacs checkpoint file in
2469      * to place and writing a FAH checkpoint is not recoverable.  Thus
2470      * the time between these operations must be kept as short a
2471      * possible.
2472      */
2473     fcCheckpoint();
2474 #endif
2475 }
2476
2477 static void check_int(FILE* fplog, const char* type, int p, int f, gmx_bool* mm)
2478 {
2479     bool foundMismatch = (p != f);
2480     if (!foundMismatch)
2481     {
2482         return;
2483     }
2484     *mm = TRUE;
2485     if (fplog)
2486     {
2487         fprintf(fplog, "  %s mismatch,\n", type);
2488         fprintf(fplog, "    current program: %d\n", p);
2489         fprintf(fplog, "    checkpoint file: %d\n", f);
2490         fprintf(fplog, "\n");
2491     }
2492 }
2493
2494 static void check_string(FILE* fplog, const char* type, const char* p, const char* f, gmx_bool* mm)
2495 {
2496     bool foundMismatch = (std::strcmp(p, f) != 0);
2497     if (!foundMismatch)
2498     {
2499         return;
2500     }
2501     *mm = TRUE;
2502     if (fplog)
2503     {
2504         fprintf(fplog, "  %s mismatch,\n", type);
2505         fprintf(fplog, "    current program: %s\n", p);
2506         fprintf(fplog, "    checkpoint file: %s\n", f);
2507         fprintf(fplog, "\n");
2508     }
2509 }
2510
2511 static void check_match(FILE*                           fplog,
2512                         const t_commrec*                cr,
2513                         const ivec                      dd_nc,
2514                         const CheckpointHeaderContents& headerContents,
2515                         gmx_bool                        reproducibilityRequested)
2516 {
2517     /* Note that this check_string on the version will also print a message
2518      * when only the minor version differs. But we only print a warning
2519      * message further down with reproducibilityRequested=TRUE.
2520      */
2521     gmx_bool versionDiffers = FALSE;
2522     check_string(fplog, "Version", gmx_version(), headerContents.version, &versionDiffers);
2523
2524     gmx_bool precisionDiffers = FALSE;
2525     check_int(fplog, "Double prec.", GMX_DOUBLE, headerContents.double_prec, &precisionDiffers);
2526     if (precisionDiffers)
2527     {
2528         const char msg_precision_difference[] =
2529                 "You are continuing a simulation with a different precision. Not matching\n"
2530                 "mixed/double precision will lead to precision or performance loss.\n";
2531         if (fplog)
2532         {
2533             fprintf(fplog, "%s\n", msg_precision_difference);
2534         }
2535     }
2536
2537     gmx_bool mm = (versionDiffers || precisionDiffers);
2538
2539     if (reproducibilityRequested)
2540     {
2541         check_string(
2542                 fplog, "Program name", gmx::getProgramContext().fullBinaryPath(), headerContents.fprog, &mm);
2543
2544         check_int(fplog, "#ranks", cr->nnodes, headerContents.nnodes, &mm);
2545     }
2546
2547     if (cr->sizeOfDefaultCommunicator > 1 && reproducibilityRequested)
2548     {
2549         // TODO: These checks are incorrect (see redmine #3309)
2550         check_int(fplog, "#PME-ranks", cr->npmenodes, headerContents.npme, &mm);
2551
2552         int npp = cr->sizeOfDefaultCommunicator;
2553         if (cr->npmenodes >= 0)
2554         {
2555             npp -= cr->npmenodes;
2556         }
2557         int npp_f = headerContents.nnodes;
2558         if (headerContents.npme >= 0)
2559         {
2560             npp_f -= headerContents.npme;
2561         }
2562         if (npp == npp_f)
2563         {
2564             check_int(fplog, "#DD-cells[x]", dd_nc[XX], headerContents.dd_nc[XX], &mm);
2565             check_int(fplog, "#DD-cells[y]", dd_nc[YY], headerContents.dd_nc[YY], &mm);
2566             check_int(fplog, "#DD-cells[z]", dd_nc[ZZ], headerContents.dd_nc[ZZ], &mm);
2567         }
2568     }
2569
2570     if (mm)
2571     {
2572         /* Gromacs should be able to continue from checkpoints between
2573          * different patch level versions, but we do not guarantee
2574          * compatibility between different major/minor versions - check this.
2575          */
2576         int gmx_major;
2577         int cpt_major;
2578         sscanf(gmx_version(), "%5d", &gmx_major);
2579         int      ret                 = sscanf(headerContents.version, "%5d", &cpt_major);
2580         gmx_bool majorVersionDiffers = (ret < 1 || gmx_major != cpt_major);
2581
2582         const char msg_major_version_difference[] =
2583                 "The current GROMACS major version is not identical to the one that\n"
2584                 "generated the checkpoint file. In principle GROMACS does not support\n"
2585                 "continuation from checkpoints between different versions, so we advise\n"
2586                 "against this. If you still want to try your luck we recommend that you use\n"
2587                 "the -noappend flag to keep your output files from the two versions separate.\n"
2588                 "This might also work around errors where the output fields in the energy\n"
2589                 "file have changed between the different versions.\n";
2590
2591         const char msg_mismatch_notice[] =
2592                 "GROMACS patchlevel, binary or parallel settings differ from previous run.\n"
2593                 "Continuation is exact, but not guaranteed to be binary identical.\n";
2594
2595         if (majorVersionDiffers)
2596         {
2597             if (fplog)
2598             {
2599                 fprintf(fplog, "%s\n", msg_major_version_difference);
2600             }
2601         }
2602         else if (reproducibilityRequested)
2603         {
2604             /* Major & minor versions match at least, but something is different. */
2605             if (fplog)
2606             {
2607                 fprintf(fplog, "%s\n", msg_mismatch_notice);
2608             }
2609         }
2610     }
2611 }
2612
2613 static void read_checkpoint(const char*                    fn,
2614                             t_fileio*                      logfio,
2615                             const t_commrec*               cr,
2616                             const ivec                     dd_nc,
2617                             IntegrationAlgorithm           eIntegrator,
2618                             int*                           init_fep_state,
2619                             CheckpointHeaderContents*      headerContents,
2620                             t_state*                       state,
2621                             ObservablesHistory*            observablesHistory,
2622                             gmx_bool                       reproducibilityRequested,
2623                             const gmx::MDModulesNotifiers& mdModulesNotifiers,
2624                             gmx::ReadCheckpointDataHolder* modularSimulatorCheckpointData,
2625                             bool                           useModularSimulator)
2626 {
2627     t_fileio* fp;
2628     char      buf[STEPSTRSIZE];
2629     int       ret;
2630
2631     fp = gmx_fio_open(fn, "r");
2632     do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, headerContents);
2633
2634     // If we are appending, then we don't want write to the open log
2635     // file because we still need to compute a checksum for it. In
2636     // that case, the filehandle will be nullptr. Otherwise, we report
2637     // to the new log file about the checkpoint file that we are
2638     // reading from.
2639     FILE* fplog = gmx_fio_getfp(logfio);
2640     if (fplog)
2641     {
2642         fprintf(fplog, "\n");
2643         fprintf(fplog, "Reading checkpoint file %s\n", fn);
2644         fprintf(fplog, "  file generated by:     %s\n", headerContents->fprog);
2645         fprintf(fplog, "  file generated at:     %s\n", headerContents->ftime);
2646         fprintf(fplog, "  GROMACS double prec.:  %d\n", headerContents->double_prec);
2647         fprintf(fplog, "  simulation part #:     %d\n", headerContents->simulation_part);
2648         fprintf(fplog, "  step:                  %s\n", gmx_step_str(headerContents->step, buf));
2649         fprintf(fplog, "  time:                  %f\n", headerContents->t);
2650         fprintf(fplog, "\n");
2651     }
2652
2653     if (headerContents->natoms != state->natoms)
2654     {
2655         gmx_fatal(FARGS,
2656                   "Checkpoint file is for a system of %d atoms, while the current system consists "
2657                   "of %d atoms",
2658                   headerContents->natoms,
2659                   state->natoms);
2660     }
2661     if (headerContents->ngtc != state->ngtc)
2662     {
2663         gmx_fatal(FARGS,
2664                   "Checkpoint file is for a system of %d T-coupling groups, while the current "
2665                   "system consists of %d T-coupling groups",
2666                   headerContents->ngtc,
2667                   state->ngtc);
2668     }
2669     if (headerContents->nnhpres != state->nnhpres)
2670     {
2671         gmx_fatal(FARGS,
2672                   "Checkpoint file is for a system of %d NH-pressure-coupling variables, while the "
2673                   "current system consists of %d NH-pressure-coupling variables",
2674                   headerContents->nnhpres,
2675                   state->nnhpres);
2676     }
2677
2678     int nlambdaHistory = (state->dfhist ? state->dfhist->nlambda : 0);
2679     if (headerContents->nlambda != nlambdaHistory)
2680     {
2681         gmx_fatal(FARGS,
2682                   "Checkpoint file is for a system with %d lambda states, while the current system "
2683                   "consists of %d lambda states",
2684                   headerContents->nlambda,
2685                   nlambdaHistory);
2686     }
2687
2688     init_gtc_state(state,
2689                    state->ngtc,
2690                    state->nnhpres,
2691                    headerContents->nhchainlength); /* need to keep this here to keep the tpr format working */
2692     /* write over whatever was read; we use the number of Nose-Hoover chains from the checkpoint */
2693
2694     if (headerContents->eIntegrator != eIntegrator)
2695     {
2696         gmx_fatal(FARGS,
2697                   "Cannot change integrator during a checkpoint restart. Perhaps you should make a "
2698                   "new .tpr with grompp -f new.mdp -t %s",
2699                   fn);
2700     }
2701
2702     // For modular simulator, no state object is populated, so we cannot do this check here!
2703     if (headerContents->flags_state != state->flags && !useModularSimulator)
2704     {
2705         gmx_fatal(FARGS,
2706                   "Cannot change a simulation algorithm during a checkpoint restart. Perhaps you "
2707                   "should make a new .tpr with grompp -f new.mdp -t %s",
2708                   fn);
2709     }
2710
2711     GMX_RELEASE_ASSERT(!(headerContents->isModularSimulatorCheckpoint && !useModularSimulator),
2712                        "Checkpoint file was written by modular simulator, but the current "
2713                        "simulation uses the legacy simulator.\n\n"
2714                        "Try the following steps:\n"
2715                        "1. Make sure the GMX_DISABLE_MODULAR_SIMULATOR environment variable is not "
2716                        "set to return to the default behavior. Retry running the simulation.\n"
2717                        "2. If the problem persists, set the environment variable "
2718                        "GMX_USE_MODULAR_SIMULATOR=ON to overwrite the default behavior and use "
2719                        "modular simulator for all implemented use cases.");
2720     GMX_RELEASE_ASSERT(!(!headerContents->isModularSimulatorCheckpoint && useModularSimulator),
2721                        "Checkpoint file was written by legacy simulator, but the current "
2722                        "simulation uses the modular simulator.\n\n"
2723                        "Try the following steps:\n"
2724                        "1. Make sure the GMX_USE_MODULAR_SIMULATOR environment variable is not set "
2725                        "to return to the default behavior. Retry running the simulation.\n"
2726                        "2. If the problem persists, set the environment variable "
2727                        "GMX_DISABLE_MODULAR_SIMULATOR=ON to overwrite the default behavior and use "
2728                        "legacy simulator for all implemented use cases.");
2729
2730     if (MASTER(cr))
2731     {
2732         check_match(fplog, cr, dd_nc, *headerContents, reproducibilityRequested);
2733     }
2734
2735     ret             = do_cpt_state(gmx_fio_getxdr(fp), headerContents->flags_state, state, nullptr);
2736     *init_fep_state = state->fep_state; /* there should be a better way to do this than setting it
2737                                            here. Investigate for 5.0. */
2738     if (ret)
2739     {
2740         cp_error();
2741     }
2742     ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents->flags_eks, &state->ekinstate, nullptr);
2743     if (ret)
2744     {
2745         cp_error();
2746     }
2747     state->ekinstate.hasReadEkinState =
2748             (((headerContents->flags_eks & enumValueToBitMask(StateKineticEntry::EkinHalfStep)) != 0)
2749              || ((headerContents->flags_eks & enumValueToBitMask(StateKineticEntry::EkinFullStep)) != 0)
2750              || ((headerContents->flags_eks & enumValueToBitMask(StateKineticEntry::EkinHalfStepOld)) != 0)
2751              || (((headerContents->flags_eks & enumValueToBitMask(StateKineticEntry::EkinNoseHooverScaleFullStep))
2752                   | (headerContents->flags_eks & enumValueToBitMask(StateKineticEntry::EkinNoseHooverScaleHalfStep))
2753                   | (headerContents->flags_eks & enumValueToBitMask(StateKineticEntry::VelocityScale)))
2754                  != 0));
2755
2756     if (headerContents->flags_enh && observablesHistory->energyHistory == nullptr)
2757     {
2758         observablesHistory->energyHistory = std::make_unique<energyhistory_t>();
2759     }
2760     ret = do_cpt_enerhist(
2761             gmx_fio_getxdr(fp), TRUE, headerContents->flags_enh, observablesHistory->energyHistory.get(), nullptr);
2762     if (ret)
2763     {
2764         cp_error();
2765     }
2766
2767     if (headerContents->flagsPullHistory)
2768     {
2769         if (observablesHistory->pullHistory == nullptr)
2770         {
2771             observablesHistory->pullHistory = std::make_unique<PullHistory>();
2772         }
2773         ret = doCptPullHist(gmx_fio_getxdr(fp),
2774                             TRUE,
2775                             headerContents->flagsPullHistory,
2776                             observablesHistory->pullHistory.get(),
2777                             nullptr);
2778         if (ret)
2779         {
2780             cp_error();
2781         }
2782     }
2783
2784     if (headerContents->file_version < CheckPointVersion::Version45)
2785     {
2786         gmx_fatal(FARGS,
2787                   "Continuing from checkpoint files written before GROMACS 4.5 is not supported");
2788     }
2789
2790     ret = do_cpt_df_hist(
2791             gmx_fio_getxdr(fp), headerContents->flags_dfh, headerContents->nlambda, &state->dfhist, nullptr);
2792     if (ret)
2793     {
2794         cp_error();
2795     }
2796
2797     if (headerContents->nED > 0 && observablesHistory->edsamHistory == nullptr)
2798     {
2799         observablesHistory->edsamHistory = std::make_unique<edsamhistory_t>(edsamhistory_t{});
2800     }
2801     ret = do_cpt_EDstate(
2802             gmx_fio_getxdr(fp), TRUE, headerContents->nED, observablesHistory->edsamHistory.get(), nullptr);
2803     if (ret)
2804     {
2805         cp_error();
2806     }
2807
2808     if (headerContents->flags_awhh != 0 && state->awhHistory == nullptr)
2809     {
2810         state->awhHistory = std::make_shared<gmx::AwhHistory>();
2811     }
2812     ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE, headerContents->flags_awhh, state->awhHistory.get(), nullptr);
2813     if (ret)
2814     {
2815         cp_error();
2816     }
2817
2818     if (headerContents->eSwapCoords != SwapType::No && observablesHistory->swapHistory == nullptr)
2819     {
2820         observablesHistory->swapHistory = std::make_unique<swaphistory_t>(swaphistory_t{});
2821     }
2822     ret = do_cpt_swapstate(
2823             gmx_fio_getxdr(fp), TRUE, headerContents->eSwapCoords, observablesHistory->swapHistory.get(), nullptr);
2824     if (ret)
2825     {
2826         cp_error();
2827     }
2828
2829     std::vector<gmx_file_position_t> outputfiles;
2830     ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, nullptr, headerContents->file_version);
2831     if (ret)
2832     {
2833         cp_error();
2834     }
2835     do_cpt_mdmodules(headerContents->file_version, fp, mdModulesNotifiers, nullptr);
2836     if (headerContents->file_version >= CheckPointVersion::ModularSimulator)
2837     {
2838         gmx::FileIOXdrSerializer serializer(fp);
2839         modularSimulatorCheckpointData->deserialize(&serializer);
2840     }
2841     ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents->file_version);
2842     if (ret)
2843     {
2844         cp_error();
2845     }
2846     if (gmx_fio_close(fp) != 0)
2847     {
2848         gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
2849     }
2850 }
2851
2852
2853 void load_checkpoint(const char*                    fn,
2854                      t_fileio*                      logfio,
2855                      const t_commrec*               cr,
2856                      const ivec                     dd_nc,
2857                      t_inputrec*                    ir,
2858                      t_state*                       state,
2859                      ObservablesHistory*            observablesHistory,
2860                      gmx_bool                       reproducibilityRequested,
2861                      const gmx::MDModulesNotifiers& mdModulesNotifiers,
2862                      gmx::ReadCheckpointDataHolder* modularSimulatorCheckpointData,
2863                      bool                           useModularSimulator)
2864 {
2865     CheckpointHeaderContents headerContents;
2866     if (SIMMASTER(cr))
2867     {
2868         /* Read the state from the checkpoint file */
2869         read_checkpoint(fn,
2870                         logfio,
2871                         cr,
2872                         dd_nc,
2873                         ir->eI,
2874                         &(ir->fepvals->init_fep_state),
2875                         &headerContents,
2876                         state,
2877                         observablesHistory,
2878                         reproducibilityRequested,
2879                         mdModulesNotifiers,
2880                         modularSimulatorCheckpointData,
2881                         useModularSimulator);
2882     }
2883     if (PAR(cr))
2884     {
2885         gmx_bcast(sizeof(headerContents.step), &headerContents.step, cr->mpiDefaultCommunicator);
2886         gmx::MDModulesCheckpointReadingBroadcast broadcastCheckPointData = {
2887             cr->mpiDefaultCommunicator, PAR(cr), headerContents.file_version
2888         };
2889         mdModulesNotifiers.checkpointingNotifier_.notify(broadcastCheckPointData);
2890     }
2891     ir->bContinuation = TRUE;
2892     if (ir->nsteps >= 0)
2893     {
2894         // TODO Should the following condition be <=? Currently if you
2895         // pass a checkpoint written by an normal completion to a restart,
2896         // mdrun will read all input, does some work but no steps, and
2897         // write successful output. But perhaps that is not desirable.
2898         // Note that we do not intend to support the use of mdrun
2899         // -nsteps to circumvent this condition.
2900         if (ir->nsteps + ir->init_step < headerContents.step)
2901         {
2902             char        buf[STEPSTRSIZE];
2903             std::string message =
2904                     gmx::formatString("The input requested %s steps, ", gmx_step_str(ir->nsteps, buf));
2905             if (ir->init_step > 0)
2906             {
2907                 message += gmx::formatString("starting from step %s, ", gmx_step_str(ir->init_step, buf));
2908             }
2909             message += gmx::formatString(
2910                     "however the checkpoint "
2911                     "file has already reached step %s. The simulation will not "
2912                     "proceed, because either your simulation is already complete, "
2913                     "or your combination of input files don't match.",
2914                     gmx_step_str(headerContents.step, buf));
2915             gmx_fatal(FARGS, "%s", message.c_str());
2916         }
2917         ir->nsteps += ir->init_step - headerContents.step;
2918     }
2919     ir->init_step       = headerContents.step;
2920     ir->simulation_part = headerContents.simulation_part + 1;
2921 }
2922
2923 void read_checkpoint_part_and_step(const char* filename, int* simulation_part, int64_t* step)
2924 {
2925     t_fileio* fp;
2926
2927     if (filename == nullptr || !gmx_fexist(filename) || ((fp = gmx_fio_open(filename, "r")) == nullptr))
2928     {
2929         *simulation_part = 0;
2930         *step            = 0;
2931         return;
2932     }
2933
2934     CheckpointHeaderContents headerContents;
2935     do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, &headerContents);
2936     gmx_fio_close(fp);
2937     *simulation_part = headerContents.simulation_part;
2938     *step            = headerContents.step;
2939 }
2940
2941 static CheckpointHeaderContents read_checkpoint_data(t_fileio*                         fp,
2942                                                      t_state*                          state,
2943                                                      std::vector<gmx_file_position_t>* outputfiles,
2944                                                      gmx::ReadCheckpointDataHolder* modularSimulatorCheckpointData)
2945 {
2946     CheckpointHeaderContents headerContents;
2947     do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, &headerContents);
2948     state->natoms        = headerContents.natoms;
2949     state->ngtc          = headerContents.ngtc;
2950     state->nnhpres       = headerContents.nnhpres;
2951     state->nhchainlength = headerContents.nhchainlength;
2952     state->flags         = headerContents.flags_state;
2953     int ret              = do_cpt_state(gmx_fio_getxdr(fp), state->flags, state, nullptr);
2954     if (ret)
2955     {
2956         cp_error();
2957     }
2958     ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents.flags_eks, &state->ekinstate, nullptr);
2959     if (ret)
2960     {
2961         cp_error();
2962     }
2963
2964     energyhistory_t enerhist;
2965     ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE, headerContents.flags_enh, &enerhist, nullptr);
2966     if (ret)
2967     {
2968         cp_error();
2969     }
2970     PullHistory pullHist = {};
2971     ret = doCptPullHist(gmx_fio_getxdr(fp), TRUE, headerContents.flagsPullHistory, &pullHist, nullptr);
2972     if (ret)
2973     {
2974         cp_error();
2975     }
2976
2977     ret = do_cpt_df_hist(
2978             gmx_fio_getxdr(fp), headerContents.flags_dfh, headerContents.nlambda, &state->dfhist, nullptr);
2979     if (ret)
2980     {
2981         cp_error();
2982     }
2983
2984     edsamhistory_t edsamhist = {};
2985     ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, headerContents.nED, &edsamhist, nullptr);
2986     if (ret)
2987     {
2988         cp_error();
2989     }
2990
2991     ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE, headerContents.flags_awhh, state->awhHistory.get(), nullptr);
2992     if (ret)
2993     {
2994         cp_error();
2995     }
2996
2997     swaphistory_t swaphist = {};
2998     ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents.eSwapCoords, &swaphist, nullptr);
2999     if (ret)
3000     {
3001         cp_error();
3002     }
3003
3004     ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, outputfiles, nullptr, headerContents.file_version);
3005
3006     if (ret)
3007     {
3008         cp_error();
3009     }
3010     gmx::MDModulesNotifiers mdModuleNotifiers;
3011     do_cpt_mdmodules(headerContents.file_version, fp, mdModuleNotifiers, nullptr);
3012     if (headerContents.file_version >= CheckPointVersion::ModularSimulator)
3013     {
3014         // Store modular checkpoint data into modularSimulatorCheckpointData
3015         gmx::FileIOXdrSerializer serializer(fp);
3016         modularSimulatorCheckpointData->deserialize(&serializer);
3017     }
3018     ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
3019     if (ret)
3020     {
3021         cp_error();
3022     }
3023     return headerContents;
3024 }
3025
3026 void read_checkpoint_trxframe(t_fileio* fp, t_trxframe* fr)
3027 {
3028     t_state                          state;
3029     std::vector<gmx_file_position_t> outputfiles;
3030     gmx::ReadCheckpointDataHolder    modularSimulatorCheckpointData;
3031     CheckpointHeaderContents         headerContents =
3032             read_checkpoint_data(fp, &state, &outputfiles, &modularSimulatorCheckpointData);
3033     if (headerContents.isModularSimulatorCheckpoint)
3034     {
3035         gmx::ModularSimulator::readCheckpointToTrxFrame(fr, &modularSimulatorCheckpointData, headerContents);
3036         return;
3037     }
3038
3039     fr->natoms    = state.natoms;
3040     fr->bStep     = TRUE;
3041     fr->step      = int64_to_int(headerContents.step, "conversion of checkpoint to trajectory");
3042     fr->bTime     = TRUE;
3043     fr->time      = headerContents.t;
3044     fr->bLambda   = TRUE;
3045     fr->lambda    = state.lambda[FreeEnergyPerturbationCouplingType::Fep];
3046     fr->fep_state = state.fep_state;
3047     fr->bAtoms    = FALSE;
3048     fr->bX        = ((state.flags & enumValueToBitMask(StateEntry::X)) != 0);
3049     if (fr->bX)
3050     {
3051         fr->x = makeRvecArray(state.x, state.natoms);
3052     }
3053     fr->bV = ((state.flags & enumValueToBitMask(StateEntry::V)) != 0);
3054     if (fr->bV)
3055     {
3056         fr->v = makeRvecArray(state.v, state.natoms);
3057     }
3058     fr->bF   = FALSE;
3059     fr->bBox = ((state.flags & enumValueToBitMask(StateEntry::Box)) != 0);
3060     if (fr->bBox)
3061     {
3062         copy_mat(state.box, fr->box);
3063     }
3064 }
3065
3066 void list_checkpoint(const char* fn, FILE* out)
3067 {
3068     t_fileio* fp;
3069     int       ret;
3070
3071     t_state state;
3072
3073     fp = gmx_fio_open(fn, "r");
3074     CheckpointHeaderContents headerContents;
3075     do_cpt_header(gmx_fio_getxdr(fp), TRUE, out, &headerContents);
3076     state.natoms        = headerContents.natoms;
3077     state.ngtc          = headerContents.ngtc;
3078     state.nnhpres       = headerContents.nnhpres;
3079     state.nhchainlength = headerContents.nhchainlength;
3080     state.flags         = headerContents.flags_state;
3081     ret                 = do_cpt_state(gmx_fio_getxdr(fp), state.flags, &state, out);
3082     if (ret)
3083     {
3084         cp_error();
3085     }
3086     ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents.flags_eks, &state.ekinstate, out);
3087     if (ret)
3088     {
3089         cp_error();
3090     }
3091
3092     energyhistory_t enerhist;
3093     ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE, headerContents.flags_enh, &enerhist, out);
3094
3095     if (ret == 0)
3096     {
3097         PullHistory pullHist = {};
3098         ret = doCptPullHist(gmx_fio_getxdr(fp), TRUE, headerContents.flagsPullHistory, &pullHist, out);
3099     }
3100
3101     if (ret == 0)
3102     {
3103         ret = do_cpt_df_hist(
3104                 gmx_fio_getxdr(fp), headerContents.flags_dfh, headerContents.nlambda, &state.dfhist, out);
3105     }
3106
3107     if (ret == 0)
3108     {
3109         edsamhistory_t edsamhist = {};
3110         ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, headerContents.nED, &edsamhist, out);
3111     }
3112
3113     if (ret == 0)
3114     {
3115         ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE, headerContents.flags_awhh, state.awhHistory.get(), out);
3116     }
3117
3118     if (ret == 0)
3119     {
3120         swaphistory_t swaphist = {};
3121         ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents.eSwapCoords, &swaphist, out);
3122     }
3123
3124     if (ret == 0)
3125     {
3126         std::vector<gmx_file_position_t> outputfiles;
3127         ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, out, headerContents.file_version);
3128     }
3129     gmx::MDModulesNotifiers mdModuleNotifiers;
3130     do_cpt_mdmodules(headerContents.file_version, fp, mdModuleNotifiers, out);
3131     if (headerContents.file_version >= CheckPointVersion::ModularSimulator)
3132     {
3133         gmx::FileIOXdrSerializer      serializer(fp);
3134         gmx::ReadCheckpointDataHolder modularSimulatorCheckpointData;
3135         modularSimulatorCheckpointData.deserialize(&serializer);
3136         modularSimulatorCheckpointData.dump(out);
3137     }
3138
3139     if (ret == 0)
3140     {
3141         ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
3142     }
3143
3144     if (ret)
3145     {
3146         cp_warning(out);
3147     }
3148     if (gmx_fio_close(fp) != 0)
3149     {
3150         gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
3151     }
3152 }
3153
3154 /* This routine cannot print tons of data, since it is called before the log file is opened. */
3155 CheckpointHeaderContents read_checkpoint_simulation_part_and_filenames(t_fileio* fp,
3156                                                                        std::vector<gmx_file_position_t>* outputfiles)
3157 {
3158     t_state                       state;
3159     gmx::ReadCheckpointDataHolder modularSimulatorCheckpointData;
3160     CheckpointHeaderContents      headerContents =
3161             read_checkpoint_data(fp, &state, outputfiles, &modularSimulatorCheckpointData);
3162     if (gmx_fio_close(fp) != 0)
3163     {
3164         gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
3165     }
3166     return headerContents;
3167 }