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