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