src/gromacs/fileio/checkpoint.cpp

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