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