2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2008,2009,2010,2011,2012,2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
5 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
6 * and including many others, as listed in the AUTHORS file in the
7 * top-level source directory and at http://www.gromacs.org.
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.
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.
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.
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.
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.
36 /* The source code in this file should be thread-safe.
37 Please keep it that way. */
40 #include "checkpoint.h"
49 #if GMX_NATIVE_WINDOWS
51 #include <sys/locking.h>
57 #include "buildinfo.h"
58 #include "gromacs/compat/make_unique.h"
59 #include "gromacs/fileio/filetypes.h"
60 #include "gromacs/fileio/gmxfio.h"
61 #include "gromacs/fileio/gmxfio-xdr.h"
62 #include "gromacs/fileio/xdr_datatype.h"
63 #include "gromacs/fileio/xdrf.h"
64 #include "gromacs/gmxlib/network.h"
65 #include "gromacs/math/vec.h"
66 #include "gromacs/math/vecdump.h"
67 #include "gromacs/math/vectypes.h"
68 #include "gromacs/mdtypes/awh-correlation-history.h"
69 #include "gromacs/mdtypes/awh-history.h"
70 #include "gromacs/mdtypes/commrec.h"
71 #include "gromacs/mdtypes/df_history.h"
72 #include "gromacs/mdtypes/edsamhistory.h"
73 #include "gromacs/mdtypes/energyhistory.h"
74 #include "gromacs/mdtypes/inputrec.h"
75 #include "gromacs/mdtypes/md_enums.h"
76 #include "gromacs/mdtypes/observableshistory.h"
77 #include "gromacs/mdtypes/state.h"
78 #include "gromacs/mdtypes/swaphistory.h"
79 #include "gromacs/trajectory/trajectoryframe.h"
80 #include "gromacs/utility/arrayref.h"
81 #include "gromacs/utility/baseversion.h"
82 #include "gromacs/utility/cstringutil.h"
83 #include "gromacs/utility/fatalerror.h"
84 #include "gromacs/utility/futil.h"
85 #include "gromacs/utility/gmxassert.h"
86 #include "gromacs/utility/int64_to_int.h"
87 #include "gromacs/utility/programcontext.h"
88 #include "gromacs/utility/smalloc.h"
89 #include "gromacs/utility/sysinfo.h"
90 #include "gromacs/utility/txtdump.h"
96 #define CPT_MAGIC1 171817
97 #define CPT_MAGIC2 171819
98 #define CPTSTRLEN 1024
100 /* cpt_version should normally only be changed
101 * when the header or footer format changes.
102 * The state data format itself is backward and forward compatible.
103 * But old code can not read a new entry that is present in the file
104 * (but can read a new format when new entries are not present).
106 static const int cpt_version = 17;
109 const char *est_names[estNR] =
112 "box", "box-rel", "box-v", "pres_prev",
113 "nosehoover-xi", "thermostat-integral",
114 "x", "v", "sdx-unsupported", "CGp", "LD-rng-unsupported", "LD-rng-i-unsupported",
115 "disre_initf", "disre_rm3tav",
116 "orire_initf", "orire_Dtav",
117 "svir_prev", "nosehoover-vxi", "v_eta", "vol0", "nhpres_xi", "nhpres_vxi", "fvir_prev", "fep_state", "MC-rng-unsupported", "MC-rng-i-unsupported",
122 eeksEKIN_N, eeksEKINH, eeksDEKINDL, eeksMVCOS, eeksEKINF, eeksEKINO, eeksEKINSCALEF, eeksEKINSCALEH, eeksVSCALE, eeksEKINTOTAL, eeksNR
125 static const char *eeks_names[eeksNR] =
127 "Ekin_n", "Ekinh", "dEkindlambda", "mv_cos",
128 "Ekinf", "Ekinh_old", "EkinScaleF_NHC", "EkinScaleH_NHC", "Vscale_NHC", "Ekin_Total"
132 eenhENERGY_N, eenhENERGY_AVER, eenhENERGY_SUM, eenhENERGY_NSUM,
133 eenhENERGY_SUM_SIM, eenhENERGY_NSUM_SIM,
134 eenhENERGY_NSTEPS, eenhENERGY_NSTEPS_SIM,
135 eenhENERGY_DELTA_H_NN,
136 eenhENERGY_DELTA_H_LIST,
137 eenhENERGY_DELTA_H_STARTTIME,
138 eenhENERGY_DELTA_H_STARTLAMBDA,
142 static const char *eenh_names[eenhNR] =
144 "energy_n", "energy_aver", "energy_sum", "energy_nsum",
145 "energy_sum_sim", "energy_nsum_sim",
146 "energy_nsteps", "energy_nsteps_sim",
148 "energy_delta_h_list",
149 "energy_delta_h_start_time",
150 "energy_delta_h_start_lambda"
153 /* free energy history variables -- need to be preserved over checkpoint */
155 edfhBEQUIL, edfhNATLAMBDA, edfhWLHISTO, edfhWLDELTA, edfhSUMWEIGHTS, edfhSUMDG, edfhSUMMINVAR, edfhSUMVAR,
156 edfhACCUMP, edfhACCUMM, edfhACCUMP2, edfhACCUMM2, edfhTIJ, edfhTIJEMP, edfhNR
158 /* free energy history variable names */
159 static const char *edfh_names[edfhNR] =
161 "bEquilibrated", "N_at_state", "Wang-Landau Histogram", "Wang-Landau Delta", "Weights", "Free Energies", "minvar", "variance",
162 "accumulated_plus", "accumulated_minus", "accumulated_plus_2", "accumulated_minus_2", "Tij", "Tij_empirical"
165 /* AWH biasing history variables */
168 eawhhEQUILIBRATEHISTOGRAM,
171 eawhhCOORDPOINT, eawhhUMBRELLAGRIDPOINT,
173 eawhhLOGSCALEDSAMPLEWEIGHT,
175 eawhhFORCECORRELATIONGRID,
179 static const char *eawhh_names[eawhhNR] =
182 "awh_equilibrateHistogram",
185 "awh_coordpoint", "awh_umbrellaGridpoint",
187 "awh_logScaledSampleWeight",
189 "awh_forceCorrelationGrid"
192 //! Higher level vector element type, only used for formatting checkpoint dumps
193 enum class CptElementType
195 integer, //!< integer
196 real, //!< float or double, not linked to precision of type real
197 real3, //!< float[3] or double[3], not linked to precision of type real
198 matrix3x3 //!< float[3][3] or double[3][3], not linked to precision of type real
201 //! \brief Parts of the checkpoint state, only used for reporting
204 microState, //!< The microstate of the simulated system
205 kineticEnergy, //!< Kinetic energy, needed for T/P-coupling state
206 energyHistory, //!< Energy observable statistics
207 freeEnergyHistory, //!< Free-energy state and observable statistics
208 accWeightHistogram //!< Accelerated weight histogram method state
211 //! \brief Return the name of a checkpoint entry based on part and part entry
212 static const char *entryName(StatePart part, int ecpt)
216 case StatePart::microState: return est_names [ecpt];
217 case StatePart::kineticEnergy: return eeks_names[ecpt];
218 case StatePart::energyHistory: return eenh_names[ecpt];
219 case StatePart::freeEnergyHistory: return edfh_names[ecpt];
220 case StatePart::accWeightHistogram: return eawhh_names[ecpt];
226 static void cp_warning(FILE *fp)
228 fprintf(fp, "\nWARNING: Checkpoint file is corrupted or truncated\n\n");
231 [[ noreturn ]] static void cp_error()
233 gmx_fatal(FARGS, "Checkpoint file corrupted/truncated, or maybe you are out of disk space?");
236 static void do_cpt_string_err(XDR *xd, const char *desc, gmx::ArrayRef<char> s, FILE *list)
238 char *data = s.data();
239 if (xdr_string(xd, &data, s.size()) == 0)
245 fprintf(list, "%s = %s\n", desc, data);
249 static int do_cpt_int(XDR *xd, const char *desc, int *i, FILE *list)
251 if (xdr_int(xd, i) == 0)
257 fprintf(list, "%s = %d\n", desc, *i);
262 static int do_cpt_u_chars(XDR *xd, const char *desc, int n, unsigned char *i, FILE *list)
266 fprintf(list, "%s = ", desc);
269 for (int j = 0; j < n && res; j++)
271 res &= xdr_u_char(xd, &i[j]);
274 fprintf(list, "%02x", i[j]);
289 static void do_cpt_int_err(XDR *xd, const char *desc, int *i, FILE *list)
291 if (do_cpt_int(xd, desc, i, list) < 0)
297 static void do_cpt_bool_err(XDR *xd, const char *desc, bool *b, FILE *list)
299 int i = static_cast<int>(*b);
301 if (do_cpt_int(xd, desc, &i, list) < 0)
309 static void do_cpt_step_err(XDR *xd, const char *desc, int64_t *i, FILE *list)
311 char buf[STEPSTRSIZE];
313 if (xdr_int64(xd, i) == 0)
319 fprintf(list, "%s = %s\n", desc, gmx_step_str(*i, buf));
323 static void do_cpt_double_err(XDR *xd, const char *desc, double *f, FILE *list)
325 if (xdr_double(xd, f) == 0)
331 fprintf(list, "%s = %f\n", desc, *f);
335 static void do_cpt_real_err(XDR *xd, real *f)
338 bool_t res = xdr_double(xd, f);
340 bool_t res = xdr_float(xd, f);
348 static void do_cpt_n_rvecs_err(XDR *xd, const char *desc, int n, rvec f[], FILE *list)
350 for (int i = 0; i < n; i++)
352 for (int j = 0; j < DIM; j++)
354 do_cpt_real_err(xd, &f[i][j]);
360 pr_rvecs(list, 0, desc, f, n);
364 template <typename T>
372 static const int value = xdr_datatype_int;
376 struct xdr_type<float>
378 static const int value = xdr_datatype_float;
382 struct xdr_type<double>
384 static const int value = xdr_datatype_double;
387 //! \brief Returns size in byte of an xdr_datatype
388 static inline unsigned int sizeOfXdrType(int xdrType)
392 case xdr_datatype_int:
394 case xdr_datatype_float:
395 return sizeof(float);
396 case xdr_datatype_double:
397 return sizeof(double);
398 default: GMX_RELEASE_ASSERT(false, "XDR data type not implemented");
404 //! \brief Returns the XDR process function for i/o of an XDR type
405 static inline xdrproc_t xdrProc(int xdrType)
409 case xdr_datatype_int:
410 return reinterpret_cast<xdrproc_t>(xdr_int);
411 case xdr_datatype_float:
412 return reinterpret_cast<xdrproc_t>(xdr_float);
413 case xdr_datatype_double:
414 return reinterpret_cast<xdrproc_t>(xdr_double);
415 default: GMX_RELEASE_ASSERT(false, "XDR data type not implemented");
421 /*! \brief Lists or only reads an xdr vector from checkpoint file
423 * When list!=NULL reads and lists the \p nf vector elements of type \p xdrType.
424 * The header for the print is set by \p part and \p ecpt.
425 * The formatting of the printing is set by \p cptElementType.
426 * When list==NULL only reads the elements.
428 static bool_t listXdrVector(XDR *xd, StatePart part, int ecpt, int nf, int xdrType,
429 FILE *list, CptElementType cptElementType)
433 const unsigned int elemSize = sizeOfXdrType(xdrType);
434 std::vector<char> data(nf*elemSize);
435 res = xdr_vector(xd, data.data(), nf, elemSize, xdrProc(xdrType));
441 case xdr_datatype_int:
442 pr_ivec(list, 0, entryName(part, ecpt), reinterpret_cast<const int *>(data.data()), nf, TRUE);
444 case xdr_datatype_float:
446 if (cptElementType == CptElementType::real3)
448 pr_rvecs(list, 0, entryName(part, ecpt), reinterpret_cast<const rvec *>(data.data()), nf/3);
453 /* Note: With double precision code dumping a single precision rvec will produce float iso rvec print, but that's a minor annoyance */
454 pr_fvec(list, 0, entryName(part, ecpt), reinterpret_cast<const float *>(data.data()), nf, TRUE);
457 case xdr_datatype_double:
459 if (cptElementType == CptElementType::real3)
461 pr_rvecs(list, 0, entryName(part, ecpt), reinterpret_cast<const rvec *>(data.data()), nf/3);
466 /* Note: With single precision code dumping a double precision rvec will produce float iso rvec print, but that's a minor annoyance */
467 pr_dvec(list, 0, entryName(part, ecpt), reinterpret_cast<const double *>(data.data()), nf, TRUE);
470 default: GMX_RELEASE_ASSERT(false, "Data type not implemented for listing");
477 //! \brief Convert a double array, typed char*, to float
478 gmx_unused static void convertArrayRealPrecision(const char *c, float *v, int n)
480 const double *d = reinterpret_cast<const double *>(c);
481 for (int i = 0; i < n; i++)
483 v[i] = static_cast<float>(d[i]);
487 //! \brief Convert a float array, typed char*, to double
488 static void convertArrayRealPrecision(const char *c, double *v, int n)
490 const float *f = reinterpret_cast<const float *>(c);
491 for (int i = 0; i < n; i++)
493 v[i] = static_cast<double>(f[i]);
497 //! \brief Generate an error for trying to convert to integer
498 static void convertArrayRealPrecision(const char gmx_unused *c, int gmx_unused *v, int gmx_unused n)
500 GMX_RELEASE_ASSERT(false, "We only expect type mismatches between float and double, not integer");
503 /*! \brief Low-level routine for reading/writing a vector of reals from/to file.
505 * This is the only routine that does the actually i/o of real vector,
506 * all other routines are intermediate level routines for specific real
507 * data types, calling this routine.
508 * Currently this routine is (too) complex, since it handles both real *
509 * and std::vector<real>. Using real * is deprecated and this routine
510 * will simplify a lot when only std::vector needs to be supported.
512 * The routine is generic to vectors with different allocators,
513 * because as part of reading a checkpoint there are vectors whose
514 * size is not known until reading has progressed far enough, so a
515 * resize method must be called.
517 * When not listing, we use either v or vector, depending on which is !=NULL.
518 * If nval >= 0, nval is used; on read this should match the passed value.
519 * If nval n<0, *nptr (with v) or vector->size() is used. On read using v,
520 * the value is stored in nptr
522 template<typename T, typename AllocatorType>
523 static int doVectorLow(XDR *xd, StatePart part, int ecpt, int sflags,
525 T **v, std::vector<T, AllocatorType> *vector,
526 FILE *list, CptElementType cptElementType)
528 GMX_RELEASE_ASSERT(list != nullptr || (v != nullptr && vector == nullptr) || (v == nullptr && vector != nullptr), "Without list, we should have exactly one of v and vector != NULL");
532 int numElemInTheFile;
537 GMX_RELEASE_ASSERT(nptr == nullptr, "With nval>=0 we should have nptr==NULL");
538 numElemInTheFile = nval;
544 GMX_RELEASE_ASSERT(nptr != nullptr, "With nval<0 we should have nptr!=NULL");
545 numElemInTheFile = *nptr;
549 numElemInTheFile = vector->size();
553 /* Read/write the vector element count */
554 res = xdr_int(xd, &numElemInTheFile);
559 /* Read/write the element data type */
560 constexpr int xdrTypeInTheCode = xdr_type<T>::value;
561 int xdrTypeInTheFile = xdrTypeInTheCode;
562 res = xdr_int(xd, &xdrTypeInTheFile);
568 if (list == nullptr && (sflags & (1 << ecpt)))
572 if (numElemInTheFile != nval)
574 gmx_fatal(FARGS, "Count mismatch for state entry %s, code count is %d, file count is %d\n", entryName(part, ecpt), nval, numElemInTheFile);
577 else if (nptr != nullptr)
579 *nptr = numElemInTheFile;
582 bool typesMatch = (xdrTypeInTheFile == xdrTypeInTheCode);
586 sprintf(buf, "mismatch for state entry %s, code precision is %s, file precision is %s",
587 entryName(part, ecpt),
588 xdr_datatype_names[xdrTypeInTheCode],
589 xdr_datatype_names[xdrTypeInTheFile]);
591 /* Matching int and real should never occur, but check anyhow */
592 if (xdrTypeInTheFile == xdr_datatype_int ||
593 xdrTypeInTheCode == xdr_datatype_int)
595 gmx_fatal(FARGS, "Type %s: incompatible checkpoint formats or corrupted checkpoint file.", buf);
604 snew(*v, numElemInTheFile);
610 /* This conditional ensures that we don't resize on write.
611 * In particular in the state where this code was written
612 * PaddedRVecVector has a size of numElemInThefile and we
613 * don't want to lose that padding here.
615 if (vector->size() < static_cast<unsigned int>(numElemInTheFile))
617 vector->resize(numElemInTheFile);
625 vChar = reinterpret_cast<char *>(vp);
629 snew(vChar, numElemInTheFile*sizeOfXdrType(xdrTypeInTheFile));
631 res = xdr_vector(xd, vChar,
632 numElemInTheFile, sizeOfXdrType(xdrTypeInTheFile),
633 xdrProc(xdrTypeInTheFile));
641 /* In the old code float-double conversion came for free.
642 * In the new code we still support it, mainly because
643 * the tip4p_continue regression test makes use of this.
644 * It's an open question if we do or don't want to allow this.
646 convertArrayRealPrecision(vChar, vp, numElemInTheFile);
652 res = listXdrVector(xd, part, ecpt, numElemInTheFile, xdrTypeInTheFile,
653 list, cptElementType);
659 //! \brief Read/Write a std::vector.
660 template <typename T>
661 static int doVector(XDR *xd, StatePart part, int ecpt, int sflags,
662 std::vector<T> *vector, FILE *list)
664 return doVectorLow<T>(xd, part, ecpt, sflags, -1, nullptr, nullptr, vector, list, CptElementType::real);
667 //! \brief Read/Write an ArrayRef<real>.
668 static int doRealArrayRef(XDR *xd, StatePart part, int ecpt, int sflags,
669 gmx::ArrayRef<real> vector, FILE *list)
671 real *v_real = vector.data();
672 return doVectorLow < real, std::allocator < real>>(xd, part, ecpt, sflags, vector.size(), nullptr, &v_real, nullptr, list, CptElementType::real);
675 //! \brief Read/Write a vector whose value_type is RVec and whose allocator is \c AllocatorType.
676 template <typename AllocatorType>
677 static int doRvecVector(XDR *xd, StatePart part, int ecpt, int sflags,
678 std::vector<gmx::RVec, AllocatorType> *v, int numAtoms, FILE *list)
680 const int numReals = numAtoms*DIM;
684 GMX_RELEASE_ASSERT(sflags & (1 << ecpt), "When not listing, the flag for the entry should be set when requesting i/o");
685 GMX_RELEASE_ASSERT(v->size() >= static_cast<size_t>(numAtoms), "v should have sufficient size for numAtoms");
687 real *v_real = (*v)[0];
689 // PaddedRVecVector is padded beyond numAtoms, we should only write
691 gmx::ArrayRef<real> ref(v_real, v_real + numReals);
693 return doRealArrayRef(xd, part, ecpt, sflags, ref, list);
697 // Use the rebind facility to change the value_type of the
698 // allocator from RVec to real.
699 using realAllocator = typename AllocatorType::template rebind<real>::other;
700 return doVectorLow<real, realAllocator>(xd, part, ecpt, sflags, numReals, nullptr, nullptr, nullptr, list, CptElementType::real);
704 /* This function stores n along with the reals for reading,
705 * but on reading it assumes that n matches the value in the checkpoint file,
706 * a fatal error is generated when this is not the case.
708 static int do_cpte_reals(XDR *xd, StatePart part, int ecpt, int sflags,
709 int n, real **v, FILE *list)
711 return doVectorLow < real, std::allocator < real>>(xd, part, ecpt, sflags, n, nullptr, v, nullptr, list, CptElementType::real);
714 /* This function does the same as do_cpte_reals,
715 * except that on reading it ignores the passed value of *n
716 * and stores the value read from the checkpoint file in *n.
718 static int do_cpte_n_reals(XDR *xd, StatePart part, int ecpt, int sflags,
719 int *n, real **v, FILE *list)
721 return doVectorLow < real, std::allocator < real>>(xd, part, ecpt, sflags, -1, n, v, nullptr, list, CptElementType::real);
724 static int do_cpte_real(XDR *xd, StatePart part, int ecpt, int sflags,
727 return doVectorLow < real, std::allocator < real>>(xd, part, ecpt, sflags, 1, nullptr, &r, nullptr, list, CptElementType::real);
730 static int do_cpte_ints(XDR *xd, StatePart part, int ecpt, int sflags,
731 int n, int **v, FILE *list)
733 return doVectorLow < int, std::allocator < int>>(xd, part, ecpt, sflags, n, nullptr, v, nullptr, list, CptElementType::integer);
736 static int do_cpte_int(XDR *xd, StatePart part, int ecpt, int sflags,
739 return do_cpte_ints(xd, part, ecpt, sflags, 1, &i, list);
742 static int do_cpte_bool(XDR *xd, StatePart part, int ecpt, int sflags,
745 int i = static_cast<int>(*b);
746 int ret = do_cpte_int(xd, part, ecpt, sflags, &i, list);
751 static int do_cpte_doubles(XDR *xd, StatePart part, int ecpt, int sflags,
752 int n, double **v, FILE *list)
754 return doVectorLow < double, std::allocator < double>>(xd, part, ecpt, sflags, n, nullptr, v, nullptr, list, CptElementType::real);
757 static int do_cpte_double(XDR *xd, StatePart part, int ecpt, int sflags,
758 double *r, FILE *list)
760 return do_cpte_doubles(xd, part, ecpt, sflags, 1, &r, list);
763 static int do_cpte_matrix(XDR *xd, StatePart part, int ecpt, int sflags,
764 matrix v, FILE *list)
770 ret = doVectorLow < real, std::allocator < real>>(xd, part, ecpt, sflags,
771 DIM*DIM, nullptr, &vr, nullptr, nullptr, CptElementType::matrix3x3);
773 if (list && ret == 0)
775 pr_rvecs(list, 0, entryName(part, ecpt), v, DIM);
782 static int do_cpte_nmatrix(XDR *xd, StatePart part, int ecpt, int sflags,
783 int n, real **v, FILE *list)
787 char name[CPTSTRLEN];
794 for (i = 0; i < n; i++)
796 reti = doVectorLow < real, std::allocator < real>>(xd, part, ecpt, sflags, n, nullptr, &(v[i]), nullptr, nullptr, CptElementType::matrix3x3);
797 if (list && reti == 0)
799 sprintf(name, "%s[%d]", entryName(part, ecpt), i);
800 pr_reals(list, 0, name, v[i], n);
810 static int do_cpte_matrices(XDR *xd, StatePart part, int ecpt, int sflags,
811 int n, matrix **v, FILE *list)
814 matrix *vp, *va = nullptr;
820 res = xdr_int(xd, &nf);
825 if (list == nullptr && nf != n)
827 gmx_fatal(FARGS, "Count mismatch for state entry %s, code count is %d, file count is %d\n", entryName(part, ecpt), n, nf);
829 if (list || !(sflags & (1<<ecpt)))
842 snew(vr, nf*DIM*DIM);
843 for (i = 0; i < nf; i++)
845 for (j = 0; j < DIM; j++)
847 for (k = 0; k < DIM; k++)
849 vr[(i*DIM+j)*DIM+k] = vp[i][j][k];
853 ret = doVectorLow < real, std::allocator < real>>(xd, part, ecpt, sflags,
854 nf*DIM*DIM, nullptr, &vr, nullptr, nullptr,
855 CptElementType::matrix3x3);
856 for (i = 0; i < nf; i++)
858 for (j = 0; j < DIM; j++)
860 for (k = 0; k < DIM; k++)
862 vp[i][j][k] = vr[(i*DIM+j)*DIM+k];
868 if (list && ret == 0)
870 for (i = 0; i < nf; i++)
872 pr_rvecs(list, 0, entryName(part, ecpt), vp[i], DIM);
883 // TODO Expand this into being a container of all data for
884 // serialization of a checkpoint, which can be stored by the caller
885 // (e.g. so that mdrun doesn't have to open the checkpoint twice).
886 // This will separate issues of allocation from those of
887 // serialization, help separate comparison from reading, and have
888 // better defined transformation functions to/from trajectory frame
890 struct CheckpointHeaderContents
893 char version[CPTSTRLEN];
894 char btime[CPTSTRLEN];
895 char buser[CPTSTRLEN];
896 char bhost[CPTSTRLEN];
898 char fprog[CPTSTRLEN];
899 char ftime[CPTSTRLEN];
921 static void do_cpt_header(XDR *xd, gmx_bool bRead, FILE *list,
922 CheckpointHeaderContents *contents)
935 res = xdr_int(xd, &magic);
938 gmx_fatal(FARGS, "The checkpoint file is empty/corrupted, or maybe you are out of disk space?");
940 if (magic != CPT_MAGIC1)
942 gmx_fatal(FARGS, "Start of file magic number mismatch, checkpoint file has %d, should be %d\n"
943 "The checkpoint file is corrupted or not a checkpoint file",
949 gmx_gethostname(fhost, 255);
951 do_cpt_string_err(xd, "GROMACS version", contents->version, list);
952 do_cpt_string_err(xd, "GROMACS build time", contents->btime, list);
953 do_cpt_string_err(xd, "GROMACS build user", contents->buser, list);
954 do_cpt_string_err(xd, "GROMACS build host", contents->bhost, list);
955 do_cpt_string_err(xd, "generating program", contents->fprog, list);
956 do_cpt_string_err(xd, "generation time", contents->ftime, list);
957 contents->file_version = cpt_version;
958 do_cpt_int_err(xd, "checkpoint file version", &contents->file_version, list);
959 if (contents->file_version > cpt_version)
961 gmx_fatal(FARGS, "Attempting to read a checkpoint file of version %d with code of version %d\n", contents->file_version, cpt_version);
963 if (contents->file_version >= 13)
965 do_cpt_int_err(xd, "GROMACS double precision", &contents->double_prec, list);
969 contents->double_prec = -1;
971 if (contents->file_version >= 12)
973 do_cpt_string_err(xd, "generating host", fhost, list);
975 do_cpt_int_err(xd, "#atoms", &contents->natoms, list);
976 do_cpt_int_err(xd, "#T-coupling groups", &contents->ngtc, list);
977 if (contents->file_version >= 10)
979 do_cpt_int_err(xd, "#Nose-Hoover T-chains", &contents->nhchainlength, list);
983 contents->nhchainlength = 1;
985 if (contents->file_version >= 11)
987 do_cpt_int_err(xd, "#Nose-Hoover T-chains for barostat ", &contents->nnhpres, list);
991 contents->nnhpres = 0;
993 if (contents->file_version >= 14)
995 do_cpt_int_err(xd, "# of total lambda states ", &contents->nlambda, list);
999 contents->nlambda = 0;
1001 do_cpt_int_err(xd, "integrator", &contents->eIntegrator, list);
1002 if (contents->file_version >= 3)
1004 do_cpt_int_err(xd, "simulation part #", &contents->simulation_part, list);
1008 contents->simulation_part = 1;
1010 if (contents->file_version >= 5)
1012 do_cpt_step_err(xd, "step", &contents->step, list);
1017 do_cpt_int_err(xd, "step", &idum, list);
1018 contents->step = static_cast<int64_t>(idum);
1020 do_cpt_double_err(xd, "t", &contents->t, list);
1021 do_cpt_int_err(xd, "#PP-ranks", &contents->nnodes, list);
1022 do_cpt_int_err(xd, "dd_nc[x]", &contents->dd_nc[XX], list);
1023 do_cpt_int_err(xd, "dd_nc[y]", &contents->dd_nc[YY], list);
1024 do_cpt_int_err(xd, "dd_nc[z]", &contents->dd_nc[ZZ], list);
1025 do_cpt_int_err(xd, "#PME-only ranks", &contents->npme, list);
1026 do_cpt_int_err(xd, "state flags", &contents->flags_state, list);
1027 if (contents->file_version >= 4)
1029 do_cpt_int_err(xd, "ekin data flags", &contents->flags_eks, list);
1030 do_cpt_int_err(xd, "energy history flags", &contents->flags_enh, list);
1034 contents->flags_eks = 0;
1035 contents->flags_enh = (contents->flags_state >> (estORIRE_DTAV+1));
1036 contents->flags_state = (contents->flags_state & ~((1<<(estORIRE_DTAV+1)) |
1037 (1<<(estORIRE_DTAV+2)) |
1038 (1<<(estORIRE_DTAV+3))));
1040 if (contents->file_version >= 14)
1042 do_cpt_int_err(xd, "df history flags", &contents->flags_dfh, list);
1046 contents->flags_dfh = 0;
1049 if (contents->file_version >= 15)
1051 do_cpt_int_err(xd, "ED data sets", &contents->nED, list);
1058 if (contents->file_version >= 16)
1060 do_cpt_int_err(xd, "swap", &contents->eSwapCoords, list);
1064 contents->eSwapCoords = eswapNO;
1067 if (contents->file_version >= 17)
1069 do_cpt_int_err(xd, "AWH history flags", &contents->flags_awhh, list);
1073 contents->flags_awhh = 0;
1077 static int do_cpt_footer(XDR *xd, int file_version)
1082 if (file_version >= 2)
1085 res = xdr_int(xd, &magic);
1090 if (magic != CPT_MAGIC2)
1099 static int do_cpt_state(XDR *xd,
1100 int fflags, t_state *state,
1104 const StatePart part = StatePart::microState;
1105 const int sflags = state->flags;
1106 // If reading, state->natoms was probably just read, so
1107 // allocations need to be managed. If writing, this won't change
1108 // anything that matters.
1109 state_change_natoms(state, state->natoms);
1110 for (int i = 0; (i < estNR && ret == 0); i++)
1112 if (fflags & (1<<i))
1116 case estLAMBDA: ret = doRealArrayRef(xd, part, i, sflags, gmx::arrayRefFromArray<real>(state->lambda.data(), state->lambda.size()), list); break;
1117 case estFEPSTATE: ret = do_cpte_int (xd, part, i, sflags, &state->fep_state, list); break;
1118 case estBOX: ret = do_cpte_matrix(xd, part, i, sflags, state->box, list); break;
1119 case estBOX_REL: ret = do_cpte_matrix(xd, part, i, sflags, state->box_rel, list); break;
1120 case estBOXV: ret = do_cpte_matrix(xd, part, i, sflags, state->boxv, list); break;
1121 case estPRES_PREV: ret = do_cpte_matrix(xd, part, i, sflags, state->pres_prev, list); break;
1122 case estSVIR_PREV: ret = do_cpte_matrix(xd, part, i, sflags, state->svir_prev, list); break;
1123 case estFVIR_PREV: ret = do_cpte_matrix(xd, part, i, sflags, state->fvir_prev, list); break;
1124 case estNH_XI: ret = doVector<double>(xd, part, i, sflags, &state->nosehoover_xi, list); break;
1125 case estNH_VXI: ret = doVector<double>(xd, part, i, sflags, &state->nosehoover_vxi, list); break;
1126 case estNHPRES_XI: ret = doVector<double>(xd, part, i, sflags, &state->nhpres_xi, list); break;
1127 case estNHPRES_VXI: ret = doVector<double>(xd, part, i, sflags, &state->nhpres_vxi, list); break;
1128 case estTHERM_INT: ret = doVector<double>(xd, part, i, sflags, &state->therm_integral, list); break;
1129 case estBAROS_INT: ret = do_cpte_double(xd, part, i, sflags, &state->baros_integral, list); break;
1130 case estVETA: ret = do_cpte_real(xd, part, i, sflags, &state->veta, list); break;
1131 case estVOL0: ret = do_cpte_real(xd, part, i, sflags, &state->vol0, list); break;
1132 case estX: ret = doRvecVector(xd, part, i, sflags, &state->x, state->natoms, list); break;
1133 case estV: ret = doRvecVector(xd, part, i, sflags, &state->v, state->natoms, list); break;
1134 /* The RNG entries are no longer written,
1135 * the next 4 lines are only for reading old files.
1137 case estLD_RNG_NOTSUPPORTED: ret = do_cpte_ints(xd, part, i, sflags, 0, nullptr, list); break;
1138 case estLD_RNGI_NOTSUPPORTED: ret = do_cpte_ints(xd, part, i, sflags, 0, nullptr, list); break;
1139 case estMC_RNG_NOTSUPPORTED: ret = do_cpte_ints(xd, part, i, sflags, 0, nullptr, list); break;
1140 case estMC_RNGI_NOTSUPPORTED: ret = do_cpte_ints(xd, part, i, sflags, 0, nullptr, list); break;
1141 case estDISRE_INITF: ret = do_cpte_real (xd, part, i, sflags, &state->hist.disre_initf, list); break;
1142 case estDISRE_RM3TAV: ret = do_cpte_n_reals(xd, part, i, sflags, &state->hist.ndisrepairs, &state->hist.disre_rm3tav, list); break;
1143 case estORIRE_INITF: ret = do_cpte_real (xd, part, i, sflags, &state->hist.orire_initf, list); break;
1144 case estORIRE_DTAV: ret = do_cpte_n_reals(xd, part, i, sflags, &state->hist.norire_Dtav, &state->hist.orire_Dtav, list); break;
1146 gmx_fatal(FARGS, "Unknown state entry %d\n"
1147 "You are reading a checkpoint file written by different code, which is not supported", i);
1155 static int do_cpt_ekinstate(XDR *xd, int fflags, ekinstate_t *ekins,
1160 const StatePart part = StatePart::kineticEnergy;
1161 for (int i = 0; (i < eeksNR && ret == 0); i++)
1163 if (fflags & (1<<i))
1168 case eeksEKIN_N: ret = do_cpte_int(xd, part, i, fflags, &ekins->ekin_n, list); break;
1169 case eeksEKINH: ret = do_cpte_matrices(xd, part, i, fflags, ekins->ekin_n, &ekins->ekinh, list); break;
1170 case eeksEKINF: ret = do_cpte_matrices(xd, part, i, fflags, ekins->ekin_n, &ekins->ekinf, list); break;
1171 case eeksEKINO: ret = do_cpte_matrices(xd, part, i, fflags, ekins->ekin_n, &ekins->ekinh_old, list); break;
1172 case eeksEKINTOTAL: ret = do_cpte_matrix(xd, part, i, fflags, ekins->ekin_total, list); break;
1173 case eeksEKINSCALEF: ret = doVector<double>(xd, part, i, fflags, &ekins->ekinscalef_nhc, list); break;
1174 case eeksVSCALE: ret = doVector<double>(xd, part, i, fflags, &ekins->vscale_nhc, list); break;
1175 case eeksEKINSCALEH: ret = doVector<double>(xd, part, i, fflags, &ekins->ekinscaleh_nhc, list); break;
1176 case eeksDEKINDL: ret = do_cpte_real(xd, part, i, fflags, &ekins->dekindl, list); break;
1177 case eeksMVCOS: ret = do_cpte_real(xd, part, i, fflags, &ekins->mvcos, list); break;
1179 gmx_fatal(FARGS, "Unknown ekin data state entry %d\n"
1180 "You are probably reading a new checkpoint file with old code", i);
1189 static int do_cpt_swapstate(XDR *xd, gmx_bool bRead,
1190 int eSwapCoords, swaphistory_t *swapstate, FILE *list)
1192 int swap_cpt_version = 2;
1194 if (eSwapCoords == eswapNO)
1199 swapstate->bFromCpt = bRead;
1200 swapstate->eSwapCoords = eSwapCoords;
1202 do_cpt_int_err(xd, "swap checkpoint version", &swap_cpt_version, list);
1203 if (bRead && swap_cpt_version < 2)
1205 gmx_fatal(FARGS, "Cannot read checkpoint files that were written with old versions"
1206 "of the ion/water position swapping protocol.\n");
1209 do_cpt_int_err(xd, "swap coupling steps", &swapstate->nAverage, list);
1211 /* When reading, init_swapcoords has not been called yet,
1212 * so we have to allocate memory first. */
1213 do_cpt_int_err(xd, "number of ion types", &swapstate->nIonTypes, list);
1216 snew(swapstate->ionType, swapstate->nIonTypes);
1219 for (int ic = 0; ic < eCompNR; ic++)
1221 for (int ii = 0; ii < swapstate->nIonTypes; ii++)
1223 swapstateIons_t *gs = &swapstate->ionType[ii];
1227 do_cpt_int_err(xd, "swap requested atoms", &gs->nMolReq[ic], list);
1231 do_cpt_int_err(xd, "swap requested atoms p", gs->nMolReq_p[ic], list);
1236 do_cpt_int_err(xd, "swap influx net", &gs->inflow_net[ic], list);
1240 do_cpt_int_err(xd, "swap influx net p", gs->inflow_net_p[ic], list);
1243 if (bRead && (nullptr == gs->nMolPast[ic]) )
1245 snew(gs->nMolPast[ic], swapstate->nAverage);
1248 for (int j = 0; j < swapstate->nAverage; j++)
1252 do_cpt_int_err(xd, "swap past atom counts", &gs->nMolPast[ic][j], list);
1256 do_cpt_int_err(xd, "swap past atom counts p", &gs->nMolPast_p[ic][j], list);
1262 /* Ion flux per channel */
1263 for (int ic = 0; ic < eChanNR; ic++)
1265 for (int ii = 0; ii < swapstate->nIonTypes; ii++)
1267 swapstateIons_t *gs = &swapstate->ionType[ii];
1271 do_cpt_int_err(xd, "channel flux A->B", &gs->fluxfromAtoB[ic], list);
1275 do_cpt_int_err(xd, "channel flux A->B p", gs->fluxfromAtoB_p[ic], list);
1280 /* Ion flux leakage */
1283 do_cpt_int_err(xd, "flux leakage", &swapstate->fluxleak, list);
1287 do_cpt_int_err(xd, "flux leakage", swapstate->fluxleak_p, list);
1291 for (int ii = 0; ii < swapstate->nIonTypes; ii++)
1293 swapstateIons_t *gs = &swapstate->ionType[ii];
1295 do_cpt_int_err(xd, "number of ions", &gs->nMol, list);
1299 snew(gs->channel_label, gs->nMol);
1300 snew(gs->comp_from, gs->nMol);
1303 do_cpt_u_chars(xd, "channel history", gs->nMol, gs->channel_label, list);
1304 do_cpt_u_chars(xd, "domain history", gs->nMol, gs->comp_from, list);
1307 /* Save the last known whole positions to checkpoint
1308 * file to be able to also make multimeric channels whole in PBC */
1309 do_cpt_int_err(xd, "Ch0 atoms", &swapstate->nat[eChan0], list);
1310 do_cpt_int_err(xd, "Ch1 atoms", &swapstate->nat[eChan1], list);
1313 snew(swapstate->xc_old_whole[eChan0], swapstate->nat[eChan0]);
1314 snew(swapstate->xc_old_whole[eChan1], swapstate->nat[eChan1]);
1315 do_cpt_n_rvecs_err(xd, "Ch0 whole x", swapstate->nat[eChan0], swapstate->xc_old_whole[eChan0], list);
1316 do_cpt_n_rvecs_err(xd, "Ch1 whole x", swapstate->nat[eChan1], swapstate->xc_old_whole[eChan1], list);
1320 do_cpt_n_rvecs_err(xd, "Ch0 whole x", swapstate->nat[eChan0], *swapstate->xc_old_whole_p[eChan0], list);
1321 do_cpt_n_rvecs_err(xd, "Ch1 whole x", swapstate->nat[eChan1], *swapstate->xc_old_whole_p[eChan1], list);
1328 static int do_cpt_enerhist(XDR *xd, gmx_bool bRead,
1329 int fflags, energyhistory_t *enerhist,
1339 GMX_RELEASE_ASSERT(enerhist != nullptr, "With energy history, we need a valid enerhist pointer");
1341 /* This is stored/read for backward compatibility */
1342 int energyHistoryNumEnergies = 0;
1345 enerhist->nsteps = 0;
1347 enerhist->nsteps_sim = 0;
1348 enerhist->nsum_sim = 0;
1350 else if (enerhist != nullptr)
1352 energyHistoryNumEnergies = enerhist->ener_sum_sim.size();
1355 delta_h_history_t *deltaH = enerhist->deltaHForeignLambdas.get();
1356 const StatePart part = StatePart::energyHistory;
1357 for (int i = 0; (i < eenhNR && ret == 0); i++)
1359 if (fflags & (1<<i))
1363 case eenhENERGY_N: ret = do_cpte_int(xd, part, i, fflags, &energyHistoryNumEnergies, list); break;
1364 case eenhENERGY_AVER: ret = doVector<double>(xd, part, i, fflags, &enerhist->ener_ave, list); break;
1365 case eenhENERGY_SUM: ret = doVector<double>(xd, part, i, fflags, &enerhist->ener_sum, list); break;
1366 case eenhENERGY_NSUM: do_cpt_step_err(xd, eenh_names[i], &enerhist->nsum, list); break;
1367 case eenhENERGY_SUM_SIM: ret = doVector<double>(xd, part, i, fflags, &enerhist->ener_sum_sim, list); break;
1368 case eenhENERGY_NSUM_SIM: do_cpt_step_err(xd, eenh_names[i], &enerhist->nsum_sim, list); break;
1369 case eenhENERGY_NSTEPS: do_cpt_step_err(xd, eenh_names[i], &enerhist->nsteps, list); break;
1370 case eenhENERGY_NSTEPS_SIM: do_cpt_step_err(xd, eenh_names[i], &enerhist->nsteps_sim, list); break;
1371 case eenhENERGY_DELTA_H_NN:
1374 if (!bRead && deltaH != nullptr)
1376 numDeltaH = deltaH->dh.size();
1378 do_cpt_int_err(xd, eenh_names[i], &numDeltaH, list);
1381 if (deltaH == nullptr)
1383 enerhist->deltaHForeignLambdas = gmx::compat::make_unique<delta_h_history_t>();
1384 deltaH = enerhist->deltaHForeignLambdas.get();
1386 deltaH->dh.resize(numDeltaH);
1387 deltaH->start_lambda_set = FALSE;
1391 case eenhENERGY_DELTA_H_LIST:
1392 for (auto dh : deltaH->dh)
1394 ret = doVector<real>(xd, part, i, fflags, &dh, list);
1397 case eenhENERGY_DELTA_H_STARTTIME:
1398 ret = do_cpte_double(xd, part, i, fflags, &(deltaH->start_time), list); break;
1399 case eenhENERGY_DELTA_H_STARTLAMBDA:
1400 ret = do_cpte_double(xd, part, i, fflags, &(deltaH->start_lambda), list); break;
1402 gmx_fatal(FARGS, "Unknown energy history entry %d\n"
1403 "You are probably reading a new checkpoint file with old code", i);
1408 if ((fflags & (1<<eenhENERGY_SUM)) && !(fflags & (1<<eenhENERGY_SUM_SIM)))
1410 /* Assume we have an old file format and copy sum to sum_sim */
1411 enerhist->ener_sum_sim = enerhist->ener_sum;
1414 if ( (fflags & (1<<eenhENERGY_NSUM)) &&
1415 !(fflags & (1<<eenhENERGY_NSTEPS)))
1417 /* Assume we have an old file format and copy nsum to nsteps */
1418 enerhist->nsteps = enerhist->nsum;
1420 if ( (fflags & (1<<eenhENERGY_NSUM_SIM)) &&
1421 !(fflags & (1<<eenhENERGY_NSTEPS_SIM)))
1423 /* Assume we have an old file format and copy nsum to nsteps */
1424 enerhist->nsteps_sim = enerhist->nsum_sim;
1430 static int do_cpt_df_hist(XDR *xd, int fflags, int nlambda, df_history_t **dfhistPtr, FILE *list)
1439 if (*dfhistPtr == nullptr)
1441 snew(*dfhistPtr, 1);
1442 (*dfhistPtr)->nlambda = nlambda;
1443 init_df_history(*dfhistPtr, nlambda);
1445 df_history_t *dfhist = *dfhistPtr;
1447 const StatePart part = StatePart::freeEnergyHistory;
1448 for (int i = 0; (i < edfhNR && ret == 0); i++)
1450 if (fflags & (1<<i))
1454 case edfhBEQUIL: ret = do_cpte_bool(xd, part, i, fflags, &dfhist->bEquil, list); break;
1455 case edfhNATLAMBDA: ret = do_cpte_ints(xd, part, i, fflags, nlambda, &dfhist->n_at_lam, list); break;
1456 case edfhWLHISTO: ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->wl_histo, list); break;
1457 case edfhWLDELTA: ret = do_cpte_real(xd, part, i, fflags, &dfhist->wl_delta, list); break;
1458 case edfhSUMWEIGHTS: ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->sum_weights, list); break;
1459 case edfhSUMDG: ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->sum_dg, list); break;
1460 case edfhSUMMINVAR: ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->sum_minvar, list); break;
1461 case edfhSUMVAR: ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->sum_variance, list); break;
1462 case edfhACCUMP: ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->accum_p, list); break;
1463 case edfhACCUMM: ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->accum_m, list); break;
1464 case edfhACCUMP2: ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->accum_p2, list); break;
1465 case edfhACCUMM2: ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->accum_m2, list); break;
1466 case edfhTIJ: ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->Tij, list); break;
1467 case edfhTIJEMP: ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->Tij_empirical, list); break;
1470 gmx_fatal(FARGS, "Unknown df history entry %d\n"
1471 "You are probably reading a new checkpoint file with old code", i);
1480 /* This function stores the last whole configuration of the reference and
1481 * average structure in the .cpt file
1483 static int do_cpt_EDstate(XDR *xd, gmx_bool bRead,
1484 int nED, edsamhistory_t *EDstate, FILE *list)
1491 EDstate->bFromCpt = bRead;
1494 /* When reading, init_edsam has not been called yet,
1495 * so we have to allocate memory first. */
1498 snew(EDstate->nref, EDstate->nED);
1499 snew(EDstate->old_sref, EDstate->nED);
1500 snew(EDstate->nav, EDstate->nED);
1501 snew(EDstate->old_sav, EDstate->nED);
1504 /* Read/write the last whole conformation of SREF and SAV for each ED dataset (usually only one) */
1505 for (int i = 0; i < EDstate->nED; i++)
1509 /* Reference structure SREF */
1510 sprintf(buf, "ED%d # of atoms in reference structure", i+1);
1511 do_cpt_int_err(xd, buf, &EDstate->nref[i], list);
1512 sprintf(buf, "ED%d x_ref", i+1);
1515 snew(EDstate->old_sref[i], EDstate->nref[i]);
1516 do_cpt_n_rvecs_err(xd, buf, EDstate->nref[i], EDstate->old_sref[i], list);
1520 do_cpt_n_rvecs_err(xd, buf, EDstate->nref[i], EDstate->old_sref_p[i], list);
1523 /* Average structure SAV */
1524 sprintf(buf, "ED%d # of atoms in average structure", i+1);
1525 do_cpt_int_err(xd, buf, &EDstate->nav[i], list);
1526 sprintf(buf, "ED%d x_av", i+1);
1529 snew(EDstate->old_sav[i], EDstate->nav[i]);
1530 do_cpt_n_rvecs_err(xd, buf, EDstate->nav[i], EDstate->old_sav[i], list);
1534 do_cpt_n_rvecs_err(xd, buf, EDstate->nav[i], EDstate->old_sav_p[i], list);
1541 static int do_cpt_correlation_grid(XDR *xd, gmx_bool bRead, gmx_unused int fflags,
1542 gmx::CorrelationGridHistory *corrGrid,
1543 FILE *list, int eawhh)
1547 do_cpt_int_err(xd, eawhh_names[eawhh], &(corrGrid->numCorrelationTensors), list);
1548 do_cpt_int_err(xd, eawhh_names[eawhh], &(corrGrid->tensorSize), list);
1549 do_cpt_int_err(xd, eawhh_names[eawhh], &(corrGrid->blockDataListSize), list);
1553 initCorrelationGridHistory(corrGrid, corrGrid->numCorrelationTensors, corrGrid->tensorSize, corrGrid->blockDataListSize);
1556 for (gmx::CorrelationBlockDataHistory &blockData : corrGrid->blockDataBuffer)
1558 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumWeight), list);
1559 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumSquareWeight), list);
1560 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumWeightX), list);
1561 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumWeightY), list);
1562 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksSquareBlockWeight), list);
1563 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksBlockSquareWeight), list);
1564 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksBlockWeightBlockWeightX), list);
1565 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksBlockWeightBlockWeightY), list);
1566 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockLength), list);
1567 do_cpt_int_err(xd, eawhh_names[eawhh], &(blockData.previousBlockIndex), list);
1568 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.correlationIntegral), list);
1574 static int do_cpt_awh_bias(XDR *xd, gmx_bool bRead,
1575 int fflags, gmx::AwhBiasHistory *biasHistory,
1580 gmx::AwhBiasStateHistory *state = &biasHistory->state;
1581 for (int i = 0; (i < eawhhNR && ret == 0); i++)
1583 if (fflags & (1<<i))
1587 case eawhhIN_INITIAL:
1588 do_cpt_bool_err(xd, eawhh_names[i], &state->in_initial, list); break;
1589 case eawhhEQUILIBRATEHISTOGRAM:
1590 do_cpt_bool_err(xd, eawhh_names[i], &state->equilibrateHistogram, list); break;
1592 do_cpt_double_err(xd, eawhh_names[i], &state->histSize, list); break;
1598 numPoints = biasHistory->pointState.size();
1600 do_cpt_int_err(xd, eawhh_names[i], &numPoints, list);
1603 biasHistory->pointState.resize(numPoints);
1607 case eawhhCOORDPOINT:
1608 for (auto &psh : biasHistory->pointState)
1610 do_cpt_double_err(xd, eawhh_names[i], &psh.target, list);
1611 do_cpt_double_err(xd, eawhh_names[i], &psh.free_energy, list);
1612 do_cpt_double_err(xd, eawhh_names[i], &psh.bias, list);
1613 do_cpt_double_err(xd, eawhh_names[i], &psh.weightsum_iteration, list);
1614 do_cpt_double_err(xd, eawhh_names[i], &psh.weightsum_covering, list);
1615 do_cpt_double_err(xd, eawhh_names[i], &psh.weightsum_tot, list);
1616 do_cpt_double_err(xd, eawhh_names[i], &psh.weightsum_ref, list);
1617 do_cpt_step_err(xd, eawhh_names[i], &psh.last_update_index, list);
1618 do_cpt_double_err(xd, eawhh_names[i], &psh.log_pmfsum, list);
1619 do_cpt_double_err(xd, eawhh_names[i], &psh.visits_iteration, list);
1620 do_cpt_double_err(xd, eawhh_names[i], &psh.visits_tot, list);
1623 case eawhhUMBRELLAGRIDPOINT:
1624 do_cpt_int_err(xd, eawhh_names[i], &(state->umbrellaGridpoint), list); break;
1625 case eawhhUPDATELIST:
1626 do_cpt_int_err(xd, eawhh_names[i], &(state->origin_index_updatelist), list);
1627 do_cpt_int_err(xd, eawhh_names[i], &(state->end_index_updatelist), list);
1629 case eawhhLOGSCALEDSAMPLEWEIGHT:
1630 do_cpt_double_err(xd, eawhh_names[i], &(state->logScaledSampleWeight), list);
1631 do_cpt_double_err(xd, eawhh_names[i], &(state->maxLogScaledSampleWeight), list);
1633 case eawhhNUMUPDATES:
1634 do_cpt_step_err(xd, eawhh_names[i], &(state->numUpdates), list);
1636 case eawhhFORCECORRELATIONGRID:
1637 ret = do_cpt_correlation_grid(xd, bRead, fflags, &biasHistory->forceCorrelationGrid, list, i);
1640 gmx_fatal(FARGS, "Unknown awh history entry %d\n", i);
1648 static int do_cpt_awh(XDR *xd, gmx_bool bRead,
1649 int fflags, gmx::AwhHistory *awhHistory,
1656 std::shared_ptr<gmx::AwhHistory> awhHistoryLocal;
1658 if (awhHistory == nullptr)
1660 GMX_RELEASE_ASSERT(bRead, "do_cpt_awh should not be called for writing without an AwhHistory");
1662 awhHistoryLocal = std::make_shared<gmx::AwhHistory>();
1663 awhHistory = awhHistoryLocal.get();
1666 /* To be able to read and write the AWH data several parameters determining the layout of the AWH structs need to be known
1667 (nbias, npoints, etc.). The best thing (?) would be to have these passed to this function. When writing to a checkpoint
1668 these parameters are available in awh_history (after calling init_awh_history). When reading from a checkpoint though, there
1669 is no initialized awh_history (it is initialized and set in this function). The AWH parameters have not always been read
1670 at the time when this function is called for reading so I don't know how to pass them as input. Here, this is solved by
1671 when writing a checkpoint, also storing parameters needed for future reading of the checkpoint.
1673 Another issue is that some variables that AWH checkpoints don't have a registered enum and string (e.g. nbias below).
1674 One difficulty is the multilevel structure of the data which would need to be represented somehow. */
1679 numBias = awhHistory->bias.size();
1681 do_cpt_int_err(xd, "awh_nbias", &numBias, list);
1685 awhHistory->bias.resize(numBias);
1687 for (auto &bias : awhHistory->bias)
1689 ret = do_cpt_awh_bias(xd, bRead, fflags, &bias, list);
1695 do_cpt_double_err(xd, "awh_potential_offset", &awhHistory->potentialOffset, list);
1700 static int do_cpt_files(XDR *xd, gmx_bool bRead,
1701 std::vector<gmx_file_position_t> *outputfiles,
1702 FILE *list, int file_version)
1705 gmx_off_t mask = 0xFFFFFFFFL;
1706 int offset_high, offset_low;
1707 std::array<char, CPTSTRLEN> buf;
1708 GMX_RELEASE_ASSERT(outputfiles, "Must have valid outputfiles");
1710 // Ensure that reading pre-allocates outputfiles, while writing
1711 // writes what is already there.
1712 int nfiles = outputfiles->size();
1713 if (do_cpt_int(xd, "number of output files", &nfiles, list) != 0)
1719 outputfiles->resize(nfiles);
1722 for (auto &outputfile : *outputfiles)
1724 /* 64-bit XDR numbers are not portable, so it is stored as separate high/low fractions */
1727 do_cpt_string_err(xd, "output filename", buf, list);
1728 std::strncpy(outputfile.filename, buf.data(), buf.size()-1);
1730 if (do_cpt_int(xd, "file_offset_high", &offset_high, list) != 0)
1734 if (do_cpt_int(xd, "file_offset_low", &offset_low, list) != 0)
1738 outputfile.offset = (static_cast<gmx_off_t>(offset_high) << 32 ) | ( static_cast<gmx_off_t>(offset_low) & mask );
1742 do_cpt_string_err(xd, "output filename", outputfile.filename, list);
1744 offset = outputfile.offset;
1752 offset_low = static_cast<int>(offset & mask);
1753 offset_high = static_cast<int>((offset >> 32) & mask);
1755 if (do_cpt_int(xd, "file_offset_high", &offset_high, list) != 0)
1759 if (do_cpt_int(xd, "file_offset_low", &offset_low, list) != 0)
1764 if (file_version >= 8)
1766 if (do_cpt_int(xd, "file_checksum_size", &outputfile.chksum_size,
1771 if (do_cpt_u_chars(xd, "file_checksum", 16, outputfile.chksum, list) != 0)
1778 outputfile.chksum_size = -1;
1785 void write_checkpoint(const char *fn, gmx_bool bNumberAndKeep,
1786 FILE *fplog, const t_commrec *cr,
1787 ivec domdecCells, int nppnodes,
1788 int eIntegrator, int simulation_part,
1789 gmx_bool bExpanded, int elamstats,
1790 int64_t step, double t,
1791 t_state *state, ObservablesHistory *observablesHistory)
1794 char *fntemp; /* the temporary checkpoint file name */
1796 char buf[1024], suffix[5+STEPSTRSIZE], sbuf[STEPSTRSIZE];
1799 if (DOMAINDECOMP(cr))
1801 npmenodes = cr->npmenodes;
1809 /* make the new temporary filename */
1810 snew(fntemp, std::strlen(fn)+5+STEPSTRSIZE);
1811 std::strcpy(fntemp, fn);
1812 fntemp[std::strlen(fn) - std::strlen(ftp2ext(fn2ftp(fn))) - 1] = '\0';
1813 sprintf(suffix, "_%s%s", "step", gmx_step_str(step, sbuf));
1814 std::strcat(fntemp, suffix);
1815 std::strcat(fntemp, fn+std::strlen(fn) - std::strlen(ftp2ext(fn2ftp(fn))) - 1);
1817 /* if we can't rename, we just overwrite the cpt file.
1818 * dangerous if interrupted.
1820 snew(fntemp, std::strlen(fn));
1821 std::strcpy(fntemp, fn);
1823 char timebuf[STRLEN];
1824 gmx_format_current_time(timebuf, STRLEN);
1828 fprintf(fplog, "Writing checkpoint, step %s at %s\n\n",
1829 gmx_step_str(step, buf), timebuf);
1832 /* Get offsets for open files */
1833 auto outputfiles = gmx_fio_get_output_file_positions();
1835 fp = gmx_fio_open(fntemp, "w");
1838 if (state->ekinstate.bUpToDate)
1841 ((1<<eeksEKIN_N) | (1<<eeksEKINH) | (1<<eeksEKINF) |
1842 (1<<eeksEKINO) | (1<<eeksEKINSCALEF) | (1<<eeksEKINSCALEH) |
1843 (1<<eeksVSCALE) | (1<<eeksDEKINDL) | (1<<eeksMVCOS));
1850 energyhistory_t *enerhist = observablesHistory->energyHistory.get();
1852 if (enerhist != nullptr && (enerhist->nsum > 0 || enerhist->nsum_sim > 0))
1854 flags_enh |= (1<<eenhENERGY_N) | (1<<eenhENERGY_NSTEPS) | (1<<eenhENERGY_NSTEPS_SIM);
1855 if (enerhist->nsum > 0)
1857 flags_enh |= ((1<<eenhENERGY_AVER) | (1<<eenhENERGY_SUM) |
1858 (1<<eenhENERGY_NSUM));
1860 if (enerhist->nsum_sim > 0)
1862 flags_enh |= ((1<<eenhENERGY_SUM_SIM) | (1<<eenhENERGY_NSUM_SIM));
1864 if (enerhist->deltaHForeignLambdas != nullptr)
1866 flags_enh |= ( (1<< eenhENERGY_DELTA_H_NN) |
1867 (1<< eenhENERGY_DELTA_H_LIST) |
1868 (1<< eenhENERGY_DELTA_H_STARTTIME) |
1869 (1<< eenhENERGY_DELTA_H_STARTLAMBDA) );
1876 flags_dfh = ((1<<edfhBEQUIL) | (1<<edfhNATLAMBDA) | (1<<edfhSUMWEIGHTS) | (1<<edfhSUMDG) |
1877 (1<<edfhTIJ) | (1<<edfhTIJEMP));
1880 flags_dfh |= ((1<<edfhWLDELTA) | (1<<edfhWLHISTO));
1882 if ((elamstats == elamstatsMINVAR) || (elamstats == elamstatsBARKER) || (elamstats == elamstatsMETROPOLIS))
1884 flags_dfh |= ((1<<edfhACCUMP) | (1<<edfhACCUMM) | (1<<edfhACCUMP2) | (1<<edfhACCUMM2)
1885 | (1<<edfhSUMMINVAR) | (1<<edfhSUMVAR));
1894 if (state->awhHistory != nullptr && !state->awhHistory->bias.empty())
1896 flags_awhh |= ( (1 << eawhhIN_INITIAL) |
1897 (1 << eawhhEQUILIBRATEHISTOGRAM) |
1898 (1 << eawhhHISTSIZE) |
1899 (1 << eawhhNPOINTS) |
1900 (1 << eawhhCOORDPOINT) |
1901 (1 << eawhhUMBRELLAGRIDPOINT) |
1902 (1 << eawhhUPDATELIST) |
1903 (1 << eawhhLOGSCALEDSAMPLEWEIGHT) |
1904 (1 << eawhhNUMUPDATES) |
1905 (1 << eawhhFORCECORRELATIONGRID));
1908 /* We can check many more things now (CPU, acceleration, etc), but
1909 * it is highly unlikely to have two separate builds with exactly
1910 * the same version, user, time, and build host!
1913 int nlambda = (state->dfhist ? state->dfhist->nlambda : 0);
1915 edsamhistory_t *edsamhist = observablesHistory->edsamHistory.get();
1916 int nED = (edsamhist ? edsamhist->nED : 0);
1918 swaphistory_t *swaphist = observablesHistory->swapHistory.get();
1919 int eSwapCoords = (swaphist ? swaphist->eSwapCoords : eswapNO);
1921 CheckpointHeaderContents headerContents =
1923 0, {0}, {0}, {0}, {0}, GMX_DOUBLE, {0}, {0},
1924 eIntegrator, simulation_part, step, t, nppnodes,
1926 state->natoms, state->ngtc, state->nnhpres,
1927 state->nhchainlength, nlambda, state->flags, flags_eks, flags_enh, flags_dfh, flags_awhh,
1930 std::strcpy(headerContents.version, gmx_version());
1931 std::strcpy(headerContents.btime, BUILD_TIME);
1932 std::strcpy(headerContents.buser, BUILD_USER);
1933 std::strcpy(headerContents.bhost, BUILD_HOST);
1934 std::strcpy(headerContents.fprog, gmx::getProgramContext().fullBinaryPath());
1935 std::strcpy(headerContents.ftime, timebuf);
1936 if (DOMAINDECOMP(cr))
1938 copy_ivec(domdecCells, headerContents.dd_nc);
1941 do_cpt_header(gmx_fio_getxdr(fp), FALSE, nullptr, &headerContents);
1943 if ((do_cpt_state(gmx_fio_getxdr(fp), state->flags, state, nullptr) < 0) ||
1944 (do_cpt_ekinstate(gmx_fio_getxdr(fp), flags_eks, &state->ekinstate, nullptr) < 0) ||
1945 (do_cpt_enerhist(gmx_fio_getxdr(fp), FALSE, flags_enh, enerhist, nullptr) < 0) ||
1946 (do_cpt_df_hist(gmx_fio_getxdr(fp), flags_dfh, nlambda, &state->dfhist, nullptr) < 0) ||
1947 (do_cpt_EDstate(gmx_fio_getxdr(fp), FALSE, nED, edsamhist, nullptr) < 0) ||
1948 (do_cpt_awh(gmx_fio_getxdr(fp), FALSE, flags_awhh, state->awhHistory.get(), nullptr) < 0) ||
1949 (do_cpt_swapstate(gmx_fio_getxdr(fp), FALSE, eSwapCoords, swaphist, nullptr) < 0) ||
1950 (do_cpt_files(gmx_fio_getxdr(fp), FALSE, &outputfiles, nullptr,
1951 headerContents.file_version) < 0))
1953 gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
1956 do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
1958 /* we really, REALLY, want to make sure to physically write the checkpoint,
1959 and all the files it depends on, out to disk. Because we've
1960 opened the checkpoint with gmx_fio_open(), it's in our list
1962 ret = gmx_fio_all_output_fsync();
1968 "Cannot fsync '%s'; maybe you are out of disk space?",
1969 gmx_fio_getname(ret));
1971 if (getenv(GMX_IGNORE_FSYNC_FAILURE_ENV) == nullptr)
1977 gmx_warning("%s", buf);
1981 if (gmx_fio_close(fp) != 0)
1983 gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
1986 /* we don't move the checkpoint if the user specified they didn't want it,
1987 or if the fsyncs failed */
1989 if (!bNumberAndKeep && !ret)
1993 /* Rename the previous checkpoint file */
1994 std::strcpy(buf, fn);
1995 buf[std::strlen(fn) - std::strlen(ftp2ext(fn2ftp(fn))) - 1] = '\0';
1996 std::strcat(buf, "_prev");
1997 std::strcat(buf, fn+std::strlen(fn) - std::strlen(ftp2ext(fn2ftp(fn))) - 1);
1999 /* we copy here so that if something goes wrong between now and
2000 * the rename below, there's always a state.cpt.
2001 * If renames are atomic (such as in POSIX systems),
2002 * this copying should be unneccesary.
2004 gmx_file_copy(fn, buf, FALSE);
2005 /* We don't really care if this fails:
2006 * there's already a new checkpoint.
2009 gmx_file_rename(fn, buf);
2012 if (gmx_file_rename(fntemp, fn) != 0)
2014 gmx_file("Cannot rename checkpoint file; maybe you are out of disk space?");
2017 #endif /* GMX_NO_RENAME */
2022 /*code for alternate checkpointing scheme. moved from top of loop over
2024 fcRequestCheckPoint();
2025 if (fcCheckPointParallel( cr->nodeid, NULL, 0) == 0)
2027 gmx_fatal( 3, __FILE__, __LINE__, "Checkpoint error on step %d\n", step );
2029 #endif /* end GMX_FAHCORE block */
2032 static void check_int(FILE *fplog, const char *type, int p, int f, gmx_bool *mm)
2034 bool foundMismatch = (p != f);
2042 fprintf(fplog, " %s mismatch,\n", type);
2043 fprintf(fplog, " current program: %d\n", p);
2044 fprintf(fplog, " checkpoint file: %d\n", f);
2045 fprintf(fplog, "\n");
2049 static void check_string(FILE *fplog, const char *type, const char *p,
2050 const char *f, gmx_bool *mm)
2052 bool foundMismatch = (std::strcmp(p, f) != 0);
2060 fprintf(fplog, " %s mismatch,\n", type);
2061 fprintf(fplog, " current program: %s\n", p);
2062 fprintf(fplog, " checkpoint file: %s\n", f);
2063 fprintf(fplog, "\n");
2067 static void check_match(FILE *fplog, const t_commrec *cr, const ivec dd_nc,
2068 const CheckpointHeaderContents &headerContents,
2069 gmx_bool reproducibilityRequested)
2071 /* Note that this check_string on the version will also print a message
2072 * when only the minor version differs. But we only print a warning
2073 * message further down with reproducibilityRequested=TRUE.
2075 gmx_bool versionDiffers = FALSE;
2076 check_string(fplog, "Version", gmx_version(), headerContents.version, &versionDiffers);
2078 gmx_bool precisionDiffers = FALSE;
2079 check_int (fplog, "Double prec.", GMX_DOUBLE, headerContents.double_prec, &precisionDiffers);
2080 if (precisionDiffers)
2082 const char msg_precision_difference[] =
2083 "You are continuing a simulation with a different precision. Not matching\n"
2084 "single/double precision will lead to precision or performance loss.\n";
2087 fprintf(fplog, "%s\n", msg_precision_difference);
2091 gmx_bool mm = (versionDiffers || precisionDiffers);
2093 if (reproducibilityRequested)
2095 check_string(fplog, "Build time", BUILD_TIME, headerContents.btime, &mm);
2096 check_string(fplog, "Build user", BUILD_USER, headerContents.buser, &mm);
2097 check_string(fplog, "Build host", BUILD_HOST, headerContents.bhost, &mm);
2098 check_string(fplog, "Program name", gmx::getProgramContext().fullBinaryPath(), headerContents.fprog, &mm);
2100 check_int (fplog, "#ranks", cr->nnodes, headerContents.nnodes, &mm);
2103 if (cr->nnodes > 1 && reproducibilityRequested)
2105 check_int (fplog, "#PME-ranks", cr->npmenodes, headerContents.npme, &mm);
2107 int npp = cr->nnodes;
2108 if (cr->npmenodes >= 0)
2110 npp -= cr->npmenodes;
2112 int npp_f = headerContents.nnodes;
2113 if (headerContents.npme >= 0)
2115 npp_f -= headerContents.npme;
2119 check_int (fplog, "#DD-cells[x]", dd_nc[XX], headerContents.dd_nc[XX], &mm);
2120 check_int (fplog, "#DD-cells[y]", dd_nc[YY], headerContents.dd_nc[YY], &mm);
2121 check_int (fplog, "#DD-cells[z]", dd_nc[ZZ], headerContents.dd_nc[ZZ], &mm);
2127 /* Gromacs should be able to continue from checkpoints between
2128 * different patch level versions, but we do not guarantee
2129 * compatibility between different major/minor versions - check this.
2133 sscanf(gmx_version(), "%5d", &gmx_major);
2134 int ret = sscanf(headerContents.version, "%5d", &cpt_major);
2135 gmx_bool majorVersionDiffers = (ret < 1 || gmx_major != cpt_major);
2137 const char msg_major_version_difference[] =
2138 "The current GROMACS major version is not identical to the one that\n"
2139 "generated the checkpoint file. In principle GROMACS does not support\n"
2140 "continuation from checkpoints between different versions, so we advise\n"
2141 "against this. If you still want to try your luck we recommend that you use\n"
2142 "the -noappend flag to keep your output files from the two versions separate.\n"
2143 "This might also work around errors where the output fields in the energy\n"
2144 "file have changed between the different versions.\n";
2146 const char msg_mismatch_notice[] =
2147 "GROMACS patchlevel, binary or parallel settings differ from previous run.\n"
2148 "Continuation is exact, but not guaranteed to be binary identical.\n";
2150 if (majorVersionDiffers)
2154 fprintf(fplog, "%s\n", msg_major_version_difference);
2157 else if (reproducibilityRequested)
2159 /* Major & minor versions match at least, but something is different. */
2162 fprintf(fplog, "%s\n", msg_mismatch_notice);
2168 static void read_checkpoint(const char *fn, FILE **pfplog,
2169 const t_commrec *cr,
2172 int *init_fep_state,
2173 CheckpointHeaderContents *headerContents,
2174 t_state *state, gmx_bool *bReadEkin,
2175 ObservablesHistory *observablesHistory,
2176 gmx_bool bAppendOutputFiles, gmx_bool bForceAppend,
2177 gmx_bool reproducibilityRequested)
2181 char buf[STEPSTRSIZE];
2183 t_fileio *chksum_file;
2184 FILE * fplog = *pfplog;
2185 unsigned char digest[16];
2186 #if !defined __native_client__ && !GMX_NATIVE_WINDOWS
2187 struct flock fl; /* don't initialize here: the struct order is OS
2191 #if !defined __native_client__ && !GMX_NATIVE_WINDOWS
2192 fl.l_type = F_WRLCK;
2193 fl.l_whence = SEEK_SET;
2199 fp = gmx_fio_open(fn, "r");
2200 do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, headerContents);
2202 if (bAppendOutputFiles &&
2203 headerContents->file_version >= 13 && headerContents->double_prec != GMX_DOUBLE)
2205 gmx_fatal(FARGS, "Output file appending requested, but the code and checkpoint file precision (single/double) don't match");
2208 /* This will not be written if we do appending, since fplog is still NULL then */
2211 fprintf(fplog, "\n");
2212 fprintf(fplog, "Reading checkpoint file %s\n", fn);
2213 fprintf(fplog, " file generated by: %s\n", headerContents->fprog);
2214 fprintf(fplog, " file generated at: %s\n", headerContents->ftime);
2215 fprintf(fplog, " GROMACS build time: %s\n", headerContents->btime);
2216 fprintf(fplog, " GROMACS build user: %s\n", headerContents->buser);
2217 fprintf(fplog, " GROMACS build host: %s\n", headerContents->bhost);
2218 fprintf(fplog, " GROMACS double prec.: %d\n", headerContents->double_prec);
2219 fprintf(fplog, " simulation part #: %d\n", headerContents->simulation_part);
2220 fprintf(fplog, " step: %s\n", gmx_step_str(headerContents->step, buf));
2221 fprintf(fplog, " time: %f\n", headerContents->t);
2222 fprintf(fplog, "\n");
2225 if (headerContents->natoms != state->natoms)
2227 gmx_fatal(FARGS, "Checkpoint file is for a system of %d atoms, while the current system consists of %d atoms", headerContents->natoms, state->natoms);
2229 if (headerContents->ngtc != state->ngtc)
2231 gmx_fatal(FARGS, "Checkpoint file is for a system of %d T-coupling groups, while the current system consists of %d T-coupling groups", headerContents->ngtc, state->ngtc);
2233 if (headerContents->nnhpres != state->nnhpres)
2235 gmx_fatal(FARGS, "Checkpoint file is for a system of %d NH-pressure-coupling variables, while the current system consists of %d NH-pressure-coupling variables", headerContents->nnhpres, state->nnhpres);
2238 int nlambdaHistory = (state->dfhist ? state->dfhist->nlambda : 0);
2239 if (headerContents->nlambda != nlambdaHistory)
2241 gmx_fatal(FARGS, "Checkpoint file is for a system with %d lambda states, while the current system consists of %d lambda states", headerContents->nlambda, nlambdaHistory);
2244 init_gtc_state(state, state->ngtc, state->nnhpres, headerContents->nhchainlength); /* need to keep this here to keep the tpr format working */
2245 /* write over whatever was read; we use the number of Nose-Hoover chains from the checkpoint */
2247 if (headerContents->eIntegrator != eIntegrator)
2249 gmx_fatal(FARGS, "Cannot change integrator during a checkpoint restart. Perhaps you should make a new .tpr with grompp -f new.mdp -t %s", fn);
2252 if (headerContents->flags_state != state->flags)
2254 gmx_fatal(FARGS, "Cannot change a simulation algorithm during a checkpoint restart. Perhaps you should make a new .tpr with grompp -f new.mdp -t %s", fn);
2259 check_match(fplog, cr, dd_nc, *headerContents,
2260 reproducibilityRequested);
2263 ret = do_cpt_state(gmx_fio_getxdr(fp), headerContents->flags_state, state, nullptr);
2264 *init_fep_state = state->fep_state; /* there should be a better way to do this than setting it here.
2265 Investigate for 5.0. */
2270 ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents->flags_eks, &state->ekinstate, nullptr);
2275 *bReadEkin = (((headerContents->flags_eks & (1<<eeksEKINH)) != 0) ||
2276 ((headerContents->flags_eks & (1<<eeksEKINF)) != 0) ||
2277 ((headerContents->flags_eks & (1<<eeksEKINO)) != 0) ||
2278 (((headerContents->flags_eks & (1<<eeksEKINSCALEF)) |
2279 (headerContents->flags_eks & (1<<eeksEKINSCALEH)) |
2280 (headerContents->flags_eks & (1<<eeksVSCALE))) != 0));
2282 if (headerContents->flags_enh && observablesHistory->energyHistory == nullptr)
2284 observablesHistory->energyHistory = gmx::compat::make_unique<energyhistory_t>();
2286 ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE,
2287 headerContents->flags_enh, observablesHistory->energyHistory.get(), nullptr);
2293 if (headerContents->file_version < 6)
2295 gmx_fatal(FARGS, "Continuing from checkpoint files written before GROMACS 4.5 is not supported");
2298 ret = do_cpt_df_hist(gmx_fio_getxdr(fp), headerContents->flags_dfh, headerContents->nlambda, &state->dfhist, nullptr);
2304 if (headerContents->nED > 0 && observablesHistory->edsamHistory == nullptr)
2306 observablesHistory->edsamHistory = gmx::compat::make_unique<edsamhistory_t>(edsamhistory_t {});
2308 ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, headerContents->nED, observablesHistory->edsamHistory.get(), nullptr);
2314 if (headerContents->flags_awhh != 0 && state->awhHistory == nullptr)
2316 state->awhHistory = std::make_shared<gmx::AwhHistory>();
2318 ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE,
2319 headerContents->flags_awhh, state->awhHistory.get(), nullptr);
2325 if (headerContents->eSwapCoords != eswapNO && observablesHistory->swapHistory == nullptr)
2327 observablesHistory->swapHistory = gmx::compat::make_unique<swaphistory_t>(swaphistory_t {});
2329 ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents->eSwapCoords, observablesHistory->swapHistory.get(), nullptr);
2335 std::vector<gmx_file_position_t> outputfiles;
2336 ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, nullptr, headerContents->file_version);
2342 ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents->file_version);
2347 if (gmx_fio_close(fp) != 0)
2349 gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
2352 /* If the user wants to append to output files,
2353 * we use the file pointer positions of the output files stored
2354 * in the checkpoint file and truncate the files such that any frames
2355 * written after the checkpoint time are removed.
2356 * All files are md5sum checked such that we can be sure that
2357 * we do not truncate other (maybe imprortant) files.
2359 if (bAppendOutputFiles)
2361 if (outputfiles.empty())
2363 gmx_fatal(FARGS, "No names of output files were recorded in the checkpoint");
2365 if (fn2ftp(outputfiles[0].filename) != efLOG)
2367 /* make sure first file is log file so that it is OK to use it for
2370 gmx_fatal(FARGS, "The first output file should always be the log "
2371 "file but instead is: %s. Cannot do appending because of this condition.", outputfiles[0].filename);
2373 bool firstFile = true;
2374 for (const auto &outputfile : outputfiles)
2376 if (outputfile.offset < 0)
2378 gmx_fatal(FARGS, "The original run wrote a file called '%s' which "
2379 "is larger than 2 GB, but mdrun did not support large file"
2380 " offsets. Can not append. Run mdrun with -noappend",
2381 outputfile.filename);
2384 chksum_file = gmx_fio_open(outputfile.filename, "a");
2387 chksum_file = gmx_fio_open(outputfile.filename, "r+");
2392 /* Note that there are systems where the lock operation
2393 * will succeed, but a second process can also lock the file.
2394 * We should probably try to detect this.
2396 #if defined __native_client__
2400 #elif GMX_NATIVE_WINDOWS
2401 if (_locking(fileno(gmx_fio_getfp(chksum_file)), _LK_NBLCK, LONG_MAX) == -1)
2403 if (fcntl(fileno(gmx_fio_getfp(chksum_file)), F_SETLK, &fl) == -1)
2406 if (errno == ENOSYS)
2410 gmx_fatal(FARGS, "File locking is not supported on this system. Use -noappend or specify -append explicitly to append anyhow.");
2416 fprintf(fplog, "\nNOTE: File locking not supported on this system, will not lock %s\n\n", outputfile.filename);
2420 else if (errno == EACCES || errno == EAGAIN)
2422 gmx_fatal(FARGS, "Failed to lock: %s. Already running "
2423 "simulation?", outputfile.filename);
2427 gmx_fatal(FARGS, "Failed to lock: %s. %s.",
2428 outputfile.filename, std::strerror(errno));
2433 /* compute md5 chksum */
2434 if (outputfile.chksum_size != -1)
2436 if (gmx_fio_get_file_md5(chksum_file, outputfile.offset,
2437 digest) != outputfile.chksum_size) /*at the end of the call the file position is at the end of the file*/
2439 gmx_fatal(FARGS, "Can't read %d bytes of '%s' to compute checksum. The file has been replaced or its contents have been modified. Cannot do appending because of this condition.",
2440 outputfile.chksum_size,
2441 outputfile.filename);
2444 /* log file needs to be seeked in case we need to truncate
2445 * (other files are truncated below)*/
2448 if (gmx_fio_seek(chksum_file, outputfile.offset))
2450 gmx_fatal(FARGS, "Seek error! Failed to truncate log-file: %s.", std::strerror(errno));
2455 /* open log file here - so that lock is never lifted
2456 * after chksum is calculated */
2459 *pfplog = gmx_fio_getfp(chksum_file);
2463 gmx_fio_close(chksum_file);
2466 /* compare md5 chksum */
2467 if (outputfile.chksum_size != -1 &&
2468 memcmp(digest, outputfile.chksum, 16) != 0)
2472 fprintf(debug, "chksum for %s: ", outputfile.filename);
2473 for (j = 0; j < 16; j++)
2475 fprintf(debug, "%02x", digest[j]);
2477 fprintf(debug, "\n");
2479 gmx_fatal(FARGS, "Checksum wrong for '%s'. The file has been replaced or its contents have been modified. Cannot do appending because of this condition.",
2480 outputfile.filename);
2485 if (!firstFile) /* log file is already seeked to correct position */
2487 #if !GMX_NATIVE_WINDOWS || !defined(GMX_FAHCORE)
2488 /* For FAHCORE, we do this elsewhere*/
2489 rc = gmx_truncate(outputfile.filename, outputfile.offset);
2492 gmx_fatal(FARGS, "Truncation of file %s failed. Cannot do appending because of this failure.", outputfile.filename);
2502 void load_checkpoint(const char *fn, FILE **fplog,
2503 const t_commrec *cr, const ivec dd_nc,
2504 t_inputrec *ir, t_state *state,
2505 gmx_bool *bReadEkin,
2506 ObservablesHistory *observablesHistory,
2507 gmx_bool bAppend, gmx_bool bForceAppend,
2508 gmx_bool reproducibilityRequested)
2510 CheckpointHeaderContents headerContents;
2513 /* Read the state from the checkpoint file */
2514 read_checkpoint(fn, fplog,
2516 ir->eI, &(ir->fepvals->init_fep_state),
2518 state, bReadEkin, observablesHistory,
2519 bAppend, bForceAppend,
2520 reproducibilityRequested);
2524 gmx_bcast(sizeof(headerContents.step), &headerContents.step, cr);
2525 gmx_bcast(sizeof(*bReadEkin), bReadEkin, cr);
2527 ir->bContinuation = TRUE;
2528 if (ir->nsteps >= 0)
2530 ir->nsteps += ir->init_step - headerContents.step;
2532 ir->init_step = headerContents.step;
2533 ir->simulation_part = headerContents.simulation_part + 1;
2536 void read_checkpoint_part_and_step(const char *filename,
2537 int *simulation_part,
2542 if (filename == nullptr ||
2543 !gmx_fexist(filename) ||
2544 (!(fp = gmx_fio_open(filename, "r"))))
2546 *simulation_part = 0;
2551 CheckpointHeaderContents headerContents;
2552 do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, &headerContents);
2554 *simulation_part = headerContents.simulation_part;
2555 *step = headerContents.step;
2558 static void read_checkpoint_data(t_fileio *fp, int *simulation_part,
2559 int64_t *step, double *t, t_state *state,
2560 std::vector<gmx_file_position_t> *outputfiles)
2564 CheckpointHeaderContents headerContents;
2565 do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, &headerContents);
2566 *simulation_part = headerContents.simulation_part;
2567 *step = headerContents.step;
2568 *t = headerContents.t;
2569 state->natoms = headerContents.natoms;
2570 state->ngtc = headerContents.ngtc;
2571 state->nnhpres = headerContents.nnhpres;
2572 state->nhchainlength = headerContents.nhchainlength;
2573 state->flags = headerContents.flags_state;
2575 do_cpt_state(gmx_fio_getxdr(fp), state->flags, state, nullptr);
2580 ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents.flags_eks, &state->ekinstate, nullptr);
2586 energyhistory_t enerhist;
2587 ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE,
2588 headerContents.flags_enh, &enerhist, nullptr);
2593 ret = do_cpt_df_hist(gmx_fio_getxdr(fp), headerContents.flags_dfh, headerContents.nlambda, &state->dfhist, nullptr);
2599 edsamhistory_t edsamhist = {};
2600 ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, headerContents.nED, &edsamhist, nullptr);
2606 ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE,
2607 headerContents.flags_awhh, state->awhHistory.get(), nullptr);
2613 swaphistory_t swaphist = {};
2614 ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents.eSwapCoords, &swaphist, nullptr);
2620 ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE,
2622 nullptr, headerContents.file_version);
2629 ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
2636 void read_checkpoint_trxframe(t_fileio *fp, t_trxframe *fr)
2639 int simulation_part;
2643 std::vector<gmx_file_position_t> outputfiles;
2644 read_checkpoint_data(fp, &simulation_part, &step, &t, &state, &outputfiles);
2646 fr->natoms = state.natoms;
2648 fr->step = int64_to_int(step,
2649 "conversion of checkpoint to trajectory");
2653 fr->lambda = state.lambda[efptFEP];
2654 fr->fep_state = state.fep_state;
2656 fr->bX = ((state.flags & (1<<estX)) != 0);
2659 fr->x = makeRvecArray(state.x, state.natoms);
2661 fr->bV = ((state.flags & (1<<estV)) != 0);
2664 fr->v = makeRvecArray(state.v, state.natoms);
2667 fr->bBox = ((state.flags & (1<<estBOX)) != 0);
2670 copy_mat(state.box, fr->box);
2674 void list_checkpoint(const char *fn, FILE *out)
2681 fp = gmx_fio_open(fn, "r");
2682 CheckpointHeaderContents headerContents;
2683 do_cpt_header(gmx_fio_getxdr(fp), TRUE, out, &headerContents);
2684 state.natoms = headerContents.natoms;
2685 state.ngtc = headerContents.ngtc;
2686 state.nnhpres = headerContents.nnhpres;
2687 state.nhchainlength = headerContents.nhchainlength;
2688 state.flags = headerContents.flags_state;
2689 ret = do_cpt_state(gmx_fio_getxdr(fp), state.flags, &state, out);
2694 ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents.flags_eks, &state.ekinstate, out);
2700 energyhistory_t enerhist;
2701 ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE,
2702 headerContents.flags_enh, &enerhist, out);
2706 ret = do_cpt_df_hist(gmx_fio_getxdr(fp),
2707 headerContents.flags_dfh, headerContents.nlambda, &state.dfhist, out);
2712 edsamhistory_t edsamhist = {};
2713 ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, headerContents.nED, &edsamhist, out);
2718 ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE,
2719 headerContents.flags_awhh, state.awhHistory.get(), out);
2724 swaphistory_t swaphist = {};
2725 ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents.eSwapCoords, &swaphist, out);
2730 std::vector<gmx_file_position_t> outputfiles;
2731 ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, out, headerContents.file_version);
2736 ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
2743 if (gmx_fio_close(fp) != 0)
2745 gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
2749 /* This routine cannot print tons of data, since it is called before the log file is opened. */
2751 read_checkpoint_simulation_part_and_filenames(t_fileio *fp,
2752 int *simulation_part,
2753 std::vector<gmx_file_position_t> *outputfiles)
2759 read_checkpoint_data(fp, simulation_part, &step, &t, &state, outputfiles);
2760 if (gmx_fio_close(fp) != 0)
2762 gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");