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