src/gromacs/fileio/checkpoint.cpp

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