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