2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2008,2009,2010,2011,2012 by the GROMACS development team.
5 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
6 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
38 /* The source code in this file should be thread-safe.
39 Please keep it that way. */
42 #include "checkpoint.h"
51 #include "buildinfo.h"
52 #include "gromacs/fileio/filetypes.h"
53 #include "gromacs/fileio/gmxfio.h"
54 #include "gromacs/fileio/gmxfio_xdr.h"
55 #include "gromacs/fileio/xdr_datatype.h"
56 #include "gromacs/fileio/xdrf.h"
57 #include "gromacs/gmxlib/network.h"
58 #include "gromacs/math/vec.h"
59 #include "gromacs/math/vecdump.h"
60 #include "gromacs/math/vectypes.h"
61 #include "gromacs/mdtypes/awh_correlation_history.h"
62 #include "gromacs/mdtypes/awh_history.h"
63 #include "gromacs/mdtypes/checkpointdata.h"
64 #include "gromacs/mdtypes/commrec.h"
65 #include "gromacs/mdtypes/df_history.h"
66 #include "gromacs/mdtypes/edsamhistory.h"
67 #include "gromacs/mdtypes/energyhistory.h"
68 #include "gromacs/mdtypes/inputrec.h"
69 #include "gromacs/mdtypes/md_enums.h"
70 #include "gromacs/mdtypes/mdrunoptions.h"
71 #include "gromacs/mdtypes/observableshistory.h"
72 #include "gromacs/mdtypes/pullhistory.h"
73 #include "gromacs/mdtypes/state.h"
74 #include "gromacs/mdtypes/swaphistory.h"
75 #include "gromacs/trajectory/trajectoryframe.h"
76 #include "gromacs/utility/arrayref.h"
77 #include "gromacs/utility/baseversion.h"
78 #include "gromacs/utility/cstringutil.h"
79 #include "gromacs/utility/fatalerror.h"
80 #include "gromacs/utility/futil.h"
81 #include "gromacs/utility/gmxassert.h"
82 #include "gromacs/utility/int64_to_int.h"
83 #include "gromacs/utility/keyvaluetree.h"
84 #include "gromacs/utility/keyvaluetreebuilder.h"
85 #include "gromacs/utility/keyvaluetreeserializer.h"
86 #include "gromacs/utility/mdmodulenotification.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"
92 #define CPT_MAGIC1 171817
93 #define CPT_MAGIC2 171819
98 template<typename ValueType>
99 void readKvtCheckpointValue(compat::not_null<ValueType*> value,
100 const std::string& name,
101 const std::string& identifier,
102 const KeyValueTreeObject& kvt)
104 const std::string key = identifier + "-" + name;
105 if (!kvt.keyExists(key))
107 std::string errorMessage = "Cannot read requested checkpoint value " + key + " .";
108 GMX_THROW(InternalError(errorMessage));
110 *value = kvt[key].cast<ValueType>();
113 template void readKvtCheckpointValue(compat::not_null<std::int64_t*> value,
114 const std::string& name,
115 const std::string& identifier,
116 const KeyValueTreeObject& kvt);
117 template void readKvtCheckpointValue(compat::not_null<real*> value,
118 const std::string& name,
119 const std::string& identifier,
120 const KeyValueTreeObject& kvt);
122 template<typename ValueType>
123 void writeKvtCheckpointValue(const ValueType& value,
124 const std::string& name,
125 const std::string& identifier,
126 KeyValueTreeObjectBuilder kvtBuilder)
128 kvtBuilder.addValue<ValueType>(identifier + "-" + name, value);
131 template void writeKvtCheckpointValue(const std::int64_t& value,
132 const std::string& name,
133 const std::string& identifier,
134 KeyValueTreeObjectBuilder kvtBuilder);
135 template void writeKvtCheckpointValue(const real& value,
136 const std::string& name,
137 const std::string& identifier,
138 KeyValueTreeObjectBuilder kvtBuilder);
143 /*! \brief Enum of values that describe the contents of a cpt file
144 * whose format matches a version number
146 * The enum helps the code be more self-documenting and ensure merges
147 * do not silently resolve when two patches make the same bump. When
148 * adding new functionality, add a new element just above cptv_Count
149 * in this enumeration, and write code below that does the right thing
150 * according to the value of file_version.
154 cptv_Unknown = 17, /**< Version before numbering scheme */
155 cptv_RemoveBuildMachineInformation, /**< remove functionality that makes mdrun builds non-reproducible */
156 cptv_ComPrevStepAsPullGroupReference, /**< Allow using COM of previous step as pull group PBC reference */
157 cptv_PullAverage, /**< Added possibility to output average pull force and position */
158 cptv_MdModules, /**< Added checkpointing for MdModules */
159 cptv_ModularSimulator, /**< Added checkpointing for modular simulator */
160 cptv_Count /**< the total number of cptv versions */
163 /*! \brief Version number of the file format written to checkpoint
164 * files by this version of the code.
166 * cpt_version should normally only be changed, via adding a new field
167 * to cptv enumeration, when the header or footer format changes.
169 * The state data format itself is backward and forward compatible.
170 * But old code can not read a new entry that is present in the file
171 * (but can read a new format when new entries are not present).
173 * The cpt_version increases whenever the file format in the main
174 * development branch changes, due to an extension of the cptv enum above.
175 * Backward compatibility for reading old run input files is maintained
176 * by checking this version number against that of the file and then using
177 * the correct code path. */
178 static const int cpt_version = cptv_Count - 1;
181 const char* est_names[estNR] = { "FE-lambda",
187 "thermostat-integral",
192 "LD-rng-unsupported",
193 "LD-rng-i-unsupported",
206 "MC-rng-unsupported",
207 "MC-rng-i-unsupported",
208 "barostat-integral" };
225 static const char* eeks_names[eeksNR] = { "Ekin_n", "Ekinh", "dEkindlambda",
226 "mv_cos", "Ekinf", "Ekinh_old",
227 "EkinScaleF_NHC", "EkinScaleH_NHC", "Vscale_NHC",
239 eenhENERGY_NSTEPS_SIM,
240 eenhENERGY_DELTA_H_NN,
241 eenhENERGY_DELTA_H_LIST,
242 eenhENERGY_DELTA_H_STARTTIME,
243 eenhENERGY_DELTA_H_STARTLAMBDA,
249 epullhPULL_NUMCOORDINATES,
250 epullhPULL_NUMGROUPS,
251 epullhPULL_NUMVALUESINXSUM,
252 epullhPULL_NUMVALUESINFSUM,
258 epullcoordh_VALUE_REF_SUM,
259 epullcoordh_VALUE_SUM,
260 epullcoordh_DR01_SUM,
261 epullcoordh_DR23_SUM,
262 epullcoordh_DR45_SUM,
263 epullcoordh_FSCAL_SUM,
264 epullcoordh_DYNAX_SUM,
274 static const char* eenh_names[eenhNR] = { "energy_n",
283 "energy_delta_h_list",
284 "energy_delta_h_start_time",
285 "energy_delta_h_start_lambda" };
287 static const char* ePullhNames[epullhNR] = { "pullhistory_numcoordinates", "pullhistory_numgroups",
288 "pullhistory_numvaluesinxsum",
289 "pullhistory_numvaluesinfsum" };
291 /* free energy history variables -- need to be preserved over checkpoint */
310 /* free energy history variable names */
311 static const char* edfh_names[edfhNR] = { "bEquilibrated",
313 "Wang-Landau Histogram",
321 "accumulated_plus_2",
322 "accumulated_minus_2",
326 /* AWH biasing history variables */
330 eawhhEQUILIBRATEHISTOGRAM,
334 eawhhUMBRELLAGRIDPOINT,
336 eawhhLOGSCALEDSAMPLEWEIGHT,
338 eawhhFORCECORRELATIONGRID,
342 static const char* eawhh_names[eawhhNR] = { "awh_in_initial", "awh_equilibrateHistogram",
343 "awh_histsize", "awh_npoints",
344 "awh_coordpoint", "awh_umbrellaGridpoint",
345 "awh_updatelist", "awh_logScaledSampleWeight",
346 "awh_numupdates", "awh_forceCorrelationGrid" };
354 static const char* epull_prev_step_com_names[epullsNR] = { "Pull groups prev step COM" };
357 //! Higher level vector element type, only used for formatting checkpoint dumps
358 enum class CptElementType
360 integer, //!< integer
361 real, //!< float or double, not linked to precision of type real
362 real3, //!< float[3] or double[3], not linked to precision of type real
363 matrix3x3 //!< float[3][3] or double[3][3], not linked to precision of type real
366 //! \brief Parts of the checkpoint state, only used for reporting
369 microState, //!< The microstate of the simulated system
370 kineticEnergy, //!< Kinetic energy, needed for T/P-coupling state
371 energyHistory, //!< Energy observable statistics
372 freeEnergyHistory, //!< Free-energy state and observable statistics
373 accWeightHistogram, //!< Accelerated weight histogram method state
374 pullState, //!< COM of previous step.
375 pullHistory //!< Pull history statistics (sums since last written output)
378 //! \brief Return the name of a checkpoint entry based on part and part entry
379 static const char* entryName(StatePart part, int ecpt)
383 case StatePart::microState: return est_names[ecpt];
384 case StatePart::kineticEnergy: return eeks_names[ecpt];
385 case StatePart::energyHistory: return eenh_names[ecpt];
386 case StatePart::freeEnergyHistory: return edfh_names[ecpt];
387 case StatePart::accWeightHistogram: return eawhh_names[ecpt];
388 case StatePart::pullState: return epull_prev_step_com_names[ecpt];
389 case StatePart::pullHistory: return ePullhNames[ecpt];
395 static void cp_warning(FILE* fp)
397 fprintf(fp, "\nWARNING: Checkpoint file is corrupted or truncated\n\n");
400 [[noreturn]] static void cp_error()
402 gmx_fatal(FARGS, "Checkpoint file corrupted/truncated, or maybe you are out of disk space?");
405 static void do_cpt_string_err(XDR* xd, const char* desc, gmx::ArrayRef<char> s, FILE* list)
407 char* data = s.data();
408 if (xdr_string(xd, &data, s.size()) == 0)
414 fprintf(list, "%s = %s\n", desc, data);
418 static int do_cpt_int(XDR* xd, const char* desc, int* i, FILE* list)
420 if (xdr_int(xd, i) == 0)
426 fprintf(list, "%s = %d\n", desc, *i);
431 static int do_cpt_u_chars(XDR* xd, const char* desc, int n, unsigned char* i, FILE* list)
435 fprintf(list, "%s = ", desc);
438 for (int j = 0; j < n && res; j++)
440 res &= xdr_u_char(xd, &i[j]);
443 fprintf(list, "%02x", i[j]);
458 static void do_cpt_int_err(XDR* xd, const char* desc, int* i, FILE* list)
460 if (do_cpt_int(xd, desc, i, list) < 0)
466 static void do_cpt_bool_err(XDR* xd, const char* desc, bool* b, FILE* list)
468 int i = static_cast<int>(*b);
470 if (do_cpt_int(xd, desc, &i, list) < 0)
478 static void do_cpt_step_err(XDR* xd, const char* desc, int64_t* i, FILE* list)
480 char buf[STEPSTRSIZE];
482 if (xdr_int64(xd, i) == 0)
488 fprintf(list, "%s = %s\n", desc, gmx_step_str(*i, buf));
492 static void do_cpt_double_err(XDR* xd, const char* desc, double* f, FILE* list)
494 if (xdr_double(xd, f) == 0)
500 fprintf(list, "%s = %f\n", desc, *f);
504 static void do_cpt_real_err(XDR* xd, real* f)
507 bool_t res = xdr_double(xd, f);
509 bool_t res = xdr_float(xd, f);
517 static void do_cpt_n_rvecs_err(XDR* xd, const char* desc, int n, rvec f[], FILE* list)
519 for (int i = 0; i < n; i++)
521 for (int j = 0; j < DIM; j++)
523 do_cpt_real_err(xd, &f[i][j]);
529 pr_rvecs(list, 0, desc, f, n);
541 static const int value = xdr_datatype_int;
545 struct xdr_type<float>
547 static const int value = xdr_datatype_float;
551 struct xdr_type<double>
553 static const int value = xdr_datatype_double;
556 //! \brief Returns size in byte of an xdr_datatype
557 static inline unsigned int sizeOfXdrType(int xdrType)
561 case xdr_datatype_int: return sizeof(int);
562 case xdr_datatype_float: return sizeof(float);
563 case xdr_datatype_double: return sizeof(double);
564 default: GMX_RELEASE_ASSERT(false, "XDR data type not implemented");
570 //! \brief Returns the XDR process function for i/o of an XDR type
571 static inline xdrproc_t xdrProc(int xdrType)
575 case xdr_datatype_int: return reinterpret_cast<xdrproc_t>(xdr_int);
576 case xdr_datatype_float: return reinterpret_cast<xdrproc_t>(xdr_float);
577 case xdr_datatype_double: return reinterpret_cast<xdrproc_t>(xdr_double);
578 default: GMX_RELEASE_ASSERT(false, "XDR data type not implemented");
584 /*! \brief Lists or only reads an xdr vector from checkpoint file
586 * When list!=NULL reads and lists the \p nf vector elements of type \p xdrType.
587 * The header for the print is set by \p part and \p ecpt.
588 * The formatting of the printing is set by \p cptElementType.
589 * When list==NULL only reads the elements.
591 static bool_t listXdrVector(XDR* xd, StatePart part, int ecpt, int nf, int xdrType, FILE* list, CptElementType cptElementType)
595 const unsigned int elemSize = sizeOfXdrType(xdrType);
596 std::vector<char> data(nf * elemSize);
597 res = xdr_vector(xd, data.data(), nf, elemSize, xdrProc(xdrType));
603 case xdr_datatype_int:
604 pr_ivec(list, 0, entryName(part, ecpt), reinterpret_cast<const int*>(data.data()), nf, TRUE);
606 case xdr_datatype_float:
608 if (cptElementType == CptElementType::real3)
610 pr_rvecs(list, 0, entryName(part, ecpt),
611 reinterpret_cast<const rvec*>(data.data()), nf / 3);
616 /* Note: With double precision code dumping a single precision rvec will produce float iso rvec print, but that's a minor annoyance */
617 pr_fvec(list, 0, entryName(part, ecpt),
618 reinterpret_cast<const float*>(data.data()), nf, TRUE);
621 case xdr_datatype_double:
623 if (cptElementType == CptElementType::real3)
625 pr_rvecs(list, 0, entryName(part, ecpt),
626 reinterpret_cast<const rvec*>(data.data()), nf / 3);
631 /* Note: With single precision code dumping a double precision rvec will produce float iso rvec print, but that's a minor annoyance */
632 pr_dvec(list, 0, entryName(part, ecpt),
633 reinterpret_cast<const double*>(data.data()), nf, TRUE);
636 default: GMX_RELEASE_ASSERT(false, "Data type not implemented for listing");
643 //! \brief Convert a double array, typed char*, to float
644 gmx_unused static void convertArrayRealPrecision(const char* c, float* v, int n)
646 const double* d = reinterpret_cast<const double*>(c);
647 for (int i = 0; i < n; i++)
649 v[i] = static_cast<float>(d[i]);
653 //! \brief Convert a float array, typed char*, to double
654 static void convertArrayRealPrecision(const char* c, double* v, int n)
656 const float* f = reinterpret_cast<const float*>(c);
657 for (int i = 0; i < n; i++)
659 v[i] = static_cast<double>(f[i]);
663 //! \brief Generate an error for trying to convert to integer
664 static void convertArrayRealPrecision(const char gmx_unused* c, int gmx_unused* v, int gmx_unused n)
666 GMX_RELEASE_ASSERT(false,
667 "We only expect type mismatches between float and double, not integer");
670 /*! \brief Low-level routine for reading/writing a vector of reals from/to file.
672 * This is the only routine that does the actually i/o of real vector,
673 * all other routines are intermediate level routines for specific real
674 * data types, calling this routine.
675 * Currently this routine is (too) complex, since it handles both real *
676 * and std::vector<real>. Using real * is deprecated and this routine
677 * will simplify a lot when only std::vector needs to be supported.
679 * The routine is generic to vectors with different allocators,
680 * because as part of reading a checkpoint there are vectors whose
681 * size is not known until reading has progressed far enough, so a
682 * resize method must be called.
684 * When not listing, we use either v or vector, depending on which is !=NULL.
685 * If nval >= 0, nval is used; on read this should match the passed value.
686 * If nval n<0, *nptr (with v) or vector->size() is used. On read using v,
687 * the value is stored in nptr
689 template<typename T, typename AllocatorType>
690 static int doVectorLow(XDR* xd,
697 std::vector<T, AllocatorType>* vector,
699 CptElementType cptElementType)
701 GMX_RELEASE_ASSERT(list != nullptr || (v != nullptr && vector == nullptr)
702 || (v == nullptr && vector != nullptr),
703 "Without list, we should have exactly one of v and vector != NULL");
707 int numElemInTheFile;
712 GMX_RELEASE_ASSERT(nptr == nullptr, "With nval>=0 we should have nptr==NULL");
713 numElemInTheFile = nval;
719 GMX_RELEASE_ASSERT(nptr != nullptr, "With nval<0 we should have nptr!=NULL");
720 numElemInTheFile = *nptr;
724 numElemInTheFile = vector->size();
728 /* Read/write the vector element count */
729 res = xdr_int(xd, &numElemInTheFile);
734 /* Read/write the element data type */
735 constexpr int xdrTypeInTheCode = xdr_type<T>::value;
736 int xdrTypeInTheFile = xdrTypeInTheCode;
737 res = xdr_int(xd, &xdrTypeInTheFile);
743 if (list == nullptr && (sflags & (1 << ecpt)))
747 if (numElemInTheFile != nval)
750 "Count mismatch for state entry %s, code count is %d, file count is %d\n",
751 entryName(part, ecpt), nval, numElemInTheFile);
754 else if (nptr != nullptr)
756 *nptr = numElemInTheFile;
759 bool typesMatch = (xdrTypeInTheFile == xdrTypeInTheCode);
763 sprintf(buf, "mismatch for state entry %s, code precision is %s, file precision is %s",
764 entryName(part, ecpt), xdr_datatype_names[xdrTypeInTheCode],
765 xdr_datatype_names[xdrTypeInTheFile]);
767 /* Matching int and real should never occur, but check anyhow */
768 if (xdrTypeInTheFile == xdr_datatype_int || xdrTypeInTheCode == xdr_datatype_int)
771 "Type %s: incompatible checkpoint formats or corrupted checkpoint file.", buf);
780 snew(*v, numElemInTheFile);
786 /* This conditional ensures that we don't resize on write.
787 * In particular in the state where this code was written
788 * vector has a size of numElemInThefile and we
789 * don't want to lose that padding here.
791 if (vector->size() < static_cast<unsigned int>(numElemInTheFile))
793 vector->resize(numElemInTheFile);
801 vChar = reinterpret_cast<char*>(vp);
805 snew(vChar, numElemInTheFile * sizeOfXdrType(xdrTypeInTheFile));
807 res = xdr_vector(xd, vChar, numElemInTheFile, sizeOfXdrType(xdrTypeInTheFile),
808 xdrProc(xdrTypeInTheFile));
816 /* In the old code float-double conversion came for free.
817 * In the new code we still support it, mainly because
818 * the tip4p_continue regression test makes use of this.
819 * It's an open question if we do or don't want to allow this.
821 convertArrayRealPrecision(vChar, vp, numElemInTheFile);
827 res = listXdrVector(xd, part, ecpt, numElemInTheFile, xdrTypeInTheFile, list, cptElementType);
833 //! \brief Read/Write a std::vector, on read checks the number of elements matches \p numElements, if specified.
836 doVector(XDR* xd, StatePart part, int ecpt, int sflags, std::vector<T>* vector, FILE* list, int numElements = -1)
838 return doVectorLow<T>(xd, part, ecpt, sflags, numElements, nullptr, nullptr, vector, list,
839 CptElementType::real);
842 //! \brief Read/Write an ArrayRef<real>.
843 static int doRealArrayRef(XDR* xd, StatePart part, int ecpt, int sflags, gmx::ArrayRef<real> vector, FILE* list)
845 real* v_real = vector.data();
846 return doVectorLow<real, std::allocator<real>>(xd, part, ecpt, sflags, vector.size(), nullptr,
847 &v_real, nullptr, list, CptElementType::real);
850 //! Convert from view of RVec to view of real.
851 static gmx::ArrayRef<real> realArrayRefFromRVecArrayRef(gmx::ArrayRef<gmx::RVec> ofRVecs)
853 return gmx::arrayRefFromArray<real>(reinterpret_cast<real*>(ofRVecs.data()), ofRVecs.size() * DIM);
856 //! \brief Read/Write a PaddedVector whose value_type is RVec.
857 template<typename PaddedVectorOfRVecType>
859 doRvecVector(XDR* xd, StatePart part, int ecpt, int sflags, PaddedVectorOfRVecType* v, int numAtoms, FILE* list)
861 const int numReals = numAtoms * DIM;
866 sflags & (1 << ecpt),
867 "When not listing, the flag for the entry should be set when requesting i/o");
868 GMX_RELEASE_ASSERT(v->size() == numAtoms, "v should have sufficient size for numAtoms");
870 return doRealArrayRef(xd, part, ecpt, sflags, realArrayRefFromRVecArrayRef(makeArrayRef(*v)), list);
874 // Use the rebind facility to change the value_type of the
875 // allocator from RVec to real.
876 using realAllocator =
877 typename std::allocator_traits<typename PaddedVectorOfRVecType::allocator_type>::template rebind_alloc<real>;
878 return doVectorLow<real, realAllocator>(xd, part, ecpt, sflags, numReals, nullptr, nullptr,
879 nullptr, list, CptElementType::real);
883 /* This function stores n along with the reals for reading,
884 * but on reading it assumes that n matches the value in the checkpoint file,
885 * a fatal error is generated when this is not the case.
887 static int do_cpte_reals(XDR* xd, StatePart part, int ecpt, int sflags, int n, real** v, FILE* list)
889 return doVectorLow<real, std::allocator<real>>(xd, part, ecpt, sflags, n, nullptr, v, nullptr,
890 list, CptElementType::real);
893 /* This function does the same as do_cpte_reals,
894 * except that on reading it ignores the passed value of *n
895 * and stores the value read from the checkpoint file in *n.
897 static int do_cpte_n_reals(XDR* xd, StatePart part, int ecpt, int sflags, int* n, real** v, FILE* list)
899 return doVectorLow<real, std::allocator<real>>(xd, part, ecpt, sflags, -1, n, v, nullptr, list,
900 CptElementType::real);
903 static int do_cpte_real(XDR* xd, StatePart part, int ecpt, int sflags, real* r, FILE* list)
905 return doVectorLow<real, std::allocator<real>>(xd, part, ecpt, sflags, 1, nullptr, &r, nullptr,
906 list, CptElementType::real);
909 static int do_cpte_ints(XDR* xd, StatePart part, int ecpt, int sflags, int n, int** v, FILE* list)
911 return doVectorLow<int, std::allocator<int>>(xd, part, ecpt, sflags, n, nullptr, v, nullptr,
912 list, CptElementType::integer);
915 static int do_cpte_int(XDR* xd, StatePart part, int ecpt, int sflags, int* i, FILE* list)
917 return do_cpte_ints(xd, part, ecpt, sflags, 1, &i, list);
920 static int do_cpte_bool(XDR* xd, StatePart part, int ecpt, int sflags, bool* b, FILE* list)
922 int i = static_cast<int>(*b);
923 int ret = do_cpte_int(xd, part, ecpt, sflags, &i, list);
928 static int do_cpte_doubles(XDR* xd, StatePart part, int ecpt, int sflags, int n, double** v, FILE* list)
930 return doVectorLow<double, std::allocator<double>>(xd, part, ecpt, sflags, n, nullptr, v,
931 nullptr, list, CptElementType::real);
934 static int do_cpte_double(XDR* xd, StatePart part, int ecpt, int sflags, double* r, FILE* list)
936 return do_cpte_doubles(xd, part, ecpt, sflags, 1, &r, list);
939 static int do_cpte_matrix(XDR* xd, StatePart part, int ecpt, int sflags, matrix v, FILE* list)
945 ret = doVectorLow<real, std::allocator<real>>(xd, part, ecpt, sflags, DIM * DIM, nullptr, &vr,
946 nullptr, nullptr, CptElementType::matrix3x3);
948 if (list && ret == 0)
950 pr_rvecs(list, 0, entryName(part, ecpt), v, DIM);
957 static int do_cpte_nmatrix(XDR* xd, StatePart part, int ecpt, int sflags, int n, real** v, FILE* list)
961 char name[CPTSTRLEN];
968 for (i = 0; i < n; i++)
970 reti = doVectorLow<real, std::allocator<real>>(xd, part, ecpt, sflags, n, nullptr, &(v[i]),
971 nullptr, nullptr, CptElementType::matrix3x3);
972 if (list && reti == 0)
974 sprintf(name, "%s[%d]", entryName(part, ecpt), i);
975 pr_reals(list, 0, name, v[i], n);
985 static int do_cpte_matrices(XDR* xd, StatePart part, int ecpt, int sflags, int n, matrix** v, FILE* list)
988 matrix *vp, *va = nullptr;
994 res = xdr_int(xd, &nf);
999 if (list == nullptr && nf != n)
1001 gmx_fatal(FARGS, "Count mismatch for state entry %s, code count is %d, file count is %d\n",
1002 entryName(part, ecpt), n, nf);
1004 if (list || !(sflags & (1 << ecpt)))
1017 snew(vr, nf * DIM * DIM);
1018 for (i = 0; i < nf; i++)
1020 for (j = 0; j < DIM; j++)
1022 for (k = 0; k < DIM; k++)
1024 vr[(i * DIM + j) * DIM + k] = vp[i][j][k];
1028 ret = doVectorLow<real, std::allocator<real>>(xd, part, ecpt, sflags, nf * DIM * DIM, nullptr,
1029 &vr, nullptr, nullptr, CptElementType::matrix3x3);
1030 for (i = 0; i < nf; i++)
1032 for (j = 0; j < DIM; j++)
1034 for (k = 0; k < DIM; k++)
1036 vp[i][j][k] = vr[(i * DIM + j) * DIM + k];
1042 if (list && ret == 0)
1044 for (i = 0; i < nf; i++)
1046 pr_rvecs(list, 0, entryName(part, ecpt), vp[i], DIM);
1057 static void do_cpt_header(XDR* xd, gmx_bool bRead, FILE* list, CheckpointHeaderContents* contents)
1070 res = xdr_int(xd, &magic);
1074 "The checkpoint file is empty/corrupted, or maybe you are out of disk space?");
1076 if (magic != CPT_MAGIC1)
1079 "Start of file magic number mismatch, checkpoint file has %d, should be %d\n"
1080 "The checkpoint file is corrupted or not a checkpoint file",
1086 gmx_gethostname(fhost, 255);
1088 do_cpt_string_err(xd, "GROMACS version", contents->version, list);
1089 // The following fields are no longer ever written with meaningful
1090 // content, but because they precede the file version, there is no
1091 // good way for new code to read the old and new formats, nor a
1092 // good way for old code to avoid giving an error while reading a
1093 // new format. So we read and write a field that no longer has a
1095 do_cpt_string_err(xd, "GROMACS build time UNUSED", contents->btime_UNUSED, list);
1096 do_cpt_string_err(xd, "GROMACS build user UNUSED", contents->buser_UNUSED, list);
1097 do_cpt_string_err(xd, "GROMACS build host UNUSED", contents->bhost_UNUSED, list);
1098 do_cpt_string_err(xd, "generating program", contents->fprog, list);
1099 do_cpt_string_err(xd, "generation time", contents->ftime, list);
1100 contents->file_version = cpt_version;
1101 do_cpt_int_err(xd, "checkpoint file version", &contents->file_version, list);
1102 if (contents->file_version > cpt_version)
1105 "Attempting to read a checkpoint file of version %d with code of version %d\n",
1106 contents->file_version, cpt_version);
1108 if (contents->file_version >= 13)
1110 do_cpt_int_err(xd, "GROMACS double precision", &contents->double_prec, list);
1114 contents->double_prec = -1;
1116 if (contents->file_version >= 12)
1118 do_cpt_string_err(xd, "generating host", fhost, list);
1120 do_cpt_int_err(xd, "#atoms", &contents->natoms, list);
1121 do_cpt_int_err(xd, "#T-coupling groups", &contents->ngtc, list);
1122 if (contents->file_version >= 10)
1124 do_cpt_int_err(xd, "#Nose-Hoover T-chains", &contents->nhchainlength, list);
1128 contents->nhchainlength = 1;
1130 if (contents->file_version >= 11)
1132 do_cpt_int_err(xd, "#Nose-Hoover T-chains for barostat ", &contents->nnhpres, list);
1136 contents->nnhpres = 0;
1138 if (contents->file_version >= 14)
1140 do_cpt_int_err(xd, "# of total lambda states ", &contents->nlambda, list);
1144 contents->nlambda = 0;
1146 do_cpt_int_err(xd, "integrator", &contents->eIntegrator, list);
1147 if (contents->file_version >= 3)
1149 do_cpt_int_err(xd, "simulation part #", &contents->simulation_part, list);
1153 contents->simulation_part = 1;
1155 if (contents->file_version >= 5)
1157 do_cpt_step_err(xd, "step", &contents->step, list);
1162 do_cpt_int_err(xd, "step", &idum, list);
1163 contents->step = static_cast<int64_t>(idum);
1165 do_cpt_double_err(xd, "t", &contents->t, list);
1166 do_cpt_int_err(xd, "#PP-ranks", &contents->nnodes, list);
1167 do_cpt_int_err(xd, "dd_nc[x]", &contents->dd_nc[XX], list);
1168 do_cpt_int_err(xd, "dd_nc[y]", &contents->dd_nc[YY], list);
1169 do_cpt_int_err(xd, "dd_nc[z]", &contents->dd_nc[ZZ], list);
1170 do_cpt_int_err(xd, "#PME-only ranks", &contents->npme, list);
1171 do_cpt_int_err(xd, "state flags", &contents->flags_state, list);
1172 if (contents->file_version >= 4)
1174 do_cpt_int_err(xd, "ekin data flags", &contents->flags_eks, list);
1175 do_cpt_int_err(xd, "energy history flags", &contents->flags_enh, list);
1179 contents->flags_eks = 0;
1180 contents->flags_enh = (contents->flags_state >> (estORIRE_DTAV + 1));
1181 contents->flags_state = (contents->flags_state
1182 & ~((1 << (estORIRE_DTAV + 1)) | (1 << (estORIRE_DTAV + 2))
1183 | (1 << (estORIRE_DTAV + 3))));
1185 if (contents->file_version >= 14)
1187 do_cpt_int_err(xd, "df history flags", &contents->flags_dfh, list);
1191 contents->flags_dfh = 0;
1194 if (contents->file_version >= 15)
1196 do_cpt_int_err(xd, "ED data sets", &contents->nED, list);
1203 if (contents->file_version >= 16)
1205 do_cpt_int_err(xd, "swap", &contents->eSwapCoords, list);
1209 contents->eSwapCoords = eswapNO;
1212 if (contents->file_version >= 17)
1214 do_cpt_int_err(xd, "AWH history flags", &contents->flags_awhh, list);
1218 contents->flags_awhh = 0;
1221 if (contents->file_version >= 18)
1223 do_cpt_int_err(xd, "pull history flags", &contents->flagsPullHistory, list);
1227 contents->flagsPullHistory = 0;
1230 if (contents->file_version >= cptv_ModularSimulator)
1232 do_cpt_bool_err(xd, "Is modular simulator checkpoint",
1233 &contents->isModularSimulatorCheckpoint, list);
1237 contents->isModularSimulatorCheckpoint = false;
1241 static int do_cpt_footer(XDR* xd, int file_version)
1246 if (file_version >= 2)
1249 res = xdr_int(xd, &magic);
1254 if (magic != CPT_MAGIC2)
1263 static int do_cpt_state(XDR* xd, int fflags, t_state* state, FILE* list)
1266 const StatePart part = StatePart::microState;
1267 const int sflags = state->flags;
1268 // If reading, state->natoms was probably just read, so
1269 // allocations need to be managed. If writing, this won't change
1270 // anything that matters.
1271 state_change_natoms(state, state->natoms);
1272 for (int i = 0; (i < estNR && ret == 0); i++)
1274 if (fflags & (1 << i))
1279 ret = doRealArrayRef(
1280 xd, part, i, sflags,
1281 gmx::arrayRefFromArray<real>(state->lambda.data(), state->lambda.size()),
1285 ret = do_cpte_int(xd, part, i, sflags, &state->fep_state, list);
1287 case estBOX: ret = do_cpte_matrix(xd, part, i, sflags, state->box, list); break;
1289 ret = do_cpte_matrix(xd, part, i, sflags, state->box_rel, list);
1291 case estBOXV: ret = do_cpte_matrix(xd, part, i, sflags, state->boxv, list); break;
1293 ret = do_cpte_matrix(xd, part, i, sflags, state->pres_prev, list);
1296 ret = do_cpte_matrix(xd, part, i, sflags, state->svir_prev, list);
1299 ret = do_cpte_matrix(xd, part, i, sflags, state->fvir_prev, list);
1302 ret = doVector<double>(xd, part, i, sflags, &state->nosehoover_xi, list);
1305 ret = doVector<double>(xd, part, i, sflags, &state->nosehoover_vxi, list);
1308 ret = doVector<double>(xd, part, i, sflags, &state->nhpres_xi, list);
1311 ret = doVector<double>(xd, part, i, sflags, &state->nhpres_vxi, list);
1314 ret = doVector<double>(xd, part, i, sflags, &state->therm_integral, list);
1317 ret = do_cpte_double(xd, part, i, sflags, &state->baros_integral, list);
1319 case estVETA: ret = do_cpte_real(xd, part, i, sflags, &state->veta, list); break;
1320 case estVOL0: ret = do_cpte_real(xd, part, i, sflags, &state->vol0, list); break;
1322 ret = doRvecVector(xd, part, i, sflags, &state->x, state->natoms, list);
1325 ret = doRvecVector(xd, part, i, sflags, &state->v, state->natoms, list);
1327 /* The RNG entries are no longer written,
1328 * the next 4 lines are only for reading old files.
1329 * It's OK that three case statements fall through.
1331 case estLD_RNG_NOTSUPPORTED:
1332 case estLD_RNGI_NOTSUPPORTED:
1333 case estMC_RNG_NOTSUPPORTED:
1334 case estMC_RNGI_NOTSUPPORTED:
1335 ret = do_cpte_ints(xd, part, i, sflags, 0, nullptr, list);
1337 case estDISRE_INITF:
1338 ret = do_cpte_real(xd, part, i, sflags, &state->hist.disre_initf, list);
1340 case estDISRE_RM3TAV:
1341 ret = do_cpte_n_reals(xd, part, i, sflags, &state->hist.ndisrepairs,
1342 &state->hist.disre_rm3tav, list);
1344 case estORIRE_INITF:
1345 ret = do_cpte_real(xd, part, i, sflags, &state->hist.orire_initf, list);
1348 ret = do_cpte_n_reals(xd, part, i, sflags, &state->hist.norire_Dtav,
1349 &state->hist.orire_Dtav, list);
1351 case estPULLCOMPREVSTEP:
1352 ret = doVector<double>(xd, part, i, sflags, &state->pull_com_prev_step, list);
1356 "Unknown state entry %d\n"
1357 "You are reading a checkpoint file written by different code, which "
1366 static int do_cpt_ekinstate(XDR* xd, int fflags, ekinstate_t* ekins, FILE* list)
1370 const StatePart part = StatePart::kineticEnergy;
1371 for (int i = 0; (i < eeksNR && ret == 0); i++)
1373 if (fflags & (1 << i))
1379 ret = do_cpte_int(xd, part, i, fflags, &ekins->ekin_n, list);
1382 ret = do_cpte_matrices(xd, part, i, fflags, ekins->ekin_n, &ekins->ekinh, list);
1385 ret = do_cpte_matrices(xd, part, i, fflags, ekins->ekin_n, &ekins->ekinf, list);
1388 ret = do_cpte_matrices(xd, part, i, fflags, ekins->ekin_n, &ekins->ekinh_old, list);
1391 ret = do_cpte_matrix(xd, part, i, fflags, ekins->ekin_total, list);
1393 case eeksEKINSCALEF:
1394 ret = doVector<double>(xd, part, i, fflags, &ekins->ekinscalef_nhc, list);
1397 ret = doVector<double>(xd, part, i, fflags, &ekins->vscale_nhc, list);
1399 case eeksEKINSCALEH:
1400 ret = doVector<double>(xd, part, i, fflags, &ekins->ekinscaleh_nhc, list);
1403 ret = do_cpte_real(xd, part, i, fflags, &ekins->dekindl, list);
1405 case eeksMVCOS: ret = do_cpte_real(xd, part, i, fflags, &ekins->mvcos, list); break;
1408 "Unknown ekin data state entry %d\n"
1409 "You are probably reading a new checkpoint file with old code",
1419 static int do_cpt_swapstate(XDR* xd, gmx_bool bRead, int eSwapCoords, swaphistory_t* swapstate, FILE* list)
1421 int swap_cpt_version = 2;
1423 if (eSwapCoords == eswapNO)
1428 swapstate->bFromCpt = bRead;
1429 swapstate->eSwapCoords = eSwapCoords;
1431 do_cpt_int_err(xd, "swap checkpoint version", &swap_cpt_version, list);
1432 if (bRead && swap_cpt_version < 2)
1435 "Cannot read checkpoint files that were written with old versions"
1436 "of the ion/water position swapping protocol.\n");
1439 do_cpt_int_err(xd, "swap coupling steps", &swapstate->nAverage, list);
1441 /* When reading, init_swapcoords has not been called yet,
1442 * so we have to allocate memory first. */
1443 do_cpt_int_err(xd, "number of ion types", &swapstate->nIonTypes, list);
1446 snew(swapstate->ionType, swapstate->nIonTypes);
1449 for (int ic = 0; ic < eCompNR; ic++)
1451 for (int ii = 0; ii < swapstate->nIonTypes; ii++)
1453 swapstateIons_t* gs = &swapstate->ionType[ii];
1457 do_cpt_int_err(xd, "swap requested atoms", &gs->nMolReq[ic], list);
1461 do_cpt_int_err(xd, "swap requested atoms p", gs->nMolReq_p[ic], list);
1466 do_cpt_int_err(xd, "swap influx net", &gs->inflow_net[ic], list);
1470 do_cpt_int_err(xd, "swap influx net p", gs->inflow_net_p[ic], list);
1473 if (bRead && (nullptr == gs->nMolPast[ic]))
1475 snew(gs->nMolPast[ic], swapstate->nAverage);
1478 for (int j = 0; j < swapstate->nAverage; j++)
1482 do_cpt_int_err(xd, "swap past atom counts", &gs->nMolPast[ic][j], list);
1486 do_cpt_int_err(xd, "swap past atom counts p", &gs->nMolPast_p[ic][j], list);
1492 /* Ion flux per channel */
1493 for (int ic = 0; ic < eChanNR; ic++)
1495 for (int ii = 0; ii < swapstate->nIonTypes; ii++)
1497 swapstateIons_t* gs = &swapstate->ionType[ii];
1501 do_cpt_int_err(xd, "channel flux A->B", &gs->fluxfromAtoB[ic], list);
1505 do_cpt_int_err(xd, "channel flux A->B p", gs->fluxfromAtoB_p[ic], list);
1510 /* Ion flux leakage */
1513 do_cpt_int_err(xd, "flux leakage", &swapstate->fluxleak, list);
1517 do_cpt_int_err(xd, "flux leakage", swapstate->fluxleak_p, list);
1521 for (int ii = 0; ii < swapstate->nIonTypes; ii++)
1523 swapstateIons_t* gs = &swapstate->ionType[ii];
1525 do_cpt_int_err(xd, "number of ions", &gs->nMol, list);
1529 snew(gs->channel_label, gs->nMol);
1530 snew(gs->comp_from, gs->nMol);
1533 do_cpt_u_chars(xd, "channel history", gs->nMol, gs->channel_label, list);
1534 do_cpt_u_chars(xd, "domain history", gs->nMol, gs->comp_from, list);
1537 /* Save the last known whole positions to checkpoint
1538 * file to be able to also make multimeric channels whole in PBC */
1539 do_cpt_int_err(xd, "Ch0 atoms", &swapstate->nat[eChan0], list);
1540 do_cpt_int_err(xd, "Ch1 atoms", &swapstate->nat[eChan1], list);
1543 snew(swapstate->xc_old_whole[eChan0], swapstate->nat[eChan0]);
1544 snew(swapstate->xc_old_whole[eChan1], swapstate->nat[eChan1]);
1545 do_cpt_n_rvecs_err(xd, "Ch0 whole x", swapstate->nat[eChan0], swapstate->xc_old_whole[eChan0], list);
1546 do_cpt_n_rvecs_err(xd, "Ch1 whole x", swapstate->nat[eChan1], swapstate->xc_old_whole[eChan1], list);
1550 do_cpt_n_rvecs_err(xd, "Ch0 whole x", swapstate->nat[eChan0],
1551 *swapstate->xc_old_whole_p[eChan0], list);
1552 do_cpt_n_rvecs_err(xd, "Ch1 whole x", swapstate->nat[eChan1],
1553 *swapstate->xc_old_whole_p[eChan1], list);
1560 static int do_cpt_enerhist(XDR* xd, gmx_bool bRead, int fflags, energyhistory_t* enerhist, FILE* list)
1569 GMX_RELEASE_ASSERT(enerhist != nullptr,
1570 "With energy history, we need a valid enerhist pointer");
1572 /* This is stored/read for backward compatibility */
1573 int energyHistoryNumEnergies = 0;
1576 enerhist->nsteps = 0;
1578 enerhist->nsteps_sim = 0;
1579 enerhist->nsum_sim = 0;
1581 else if (enerhist != nullptr)
1583 energyHistoryNumEnergies = enerhist->ener_sum_sim.size();
1586 delta_h_history_t* deltaH = enerhist->deltaHForeignLambdas.get();
1587 const StatePart part = StatePart::energyHistory;
1588 for (int i = 0; (i < eenhNR && ret == 0); i++)
1590 if (fflags & (1 << i))
1595 ret = do_cpte_int(xd, part, i, fflags, &energyHistoryNumEnergies, list);
1597 case eenhENERGY_AVER:
1598 ret = doVector<double>(xd, part, i, fflags, &enerhist->ener_ave, list);
1600 case eenhENERGY_SUM:
1601 ret = doVector<double>(xd, part, i, fflags, &enerhist->ener_sum, list);
1603 case eenhENERGY_NSUM:
1604 do_cpt_step_err(xd, eenh_names[i], &enerhist->nsum, list);
1606 case eenhENERGY_SUM_SIM:
1607 ret = doVector<double>(xd, part, i, fflags, &enerhist->ener_sum_sim, list);
1609 case eenhENERGY_NSUM_SIM:
1610 do_cpt_step_err(xd, eenh_names[i], &enerhist->nsum_sim, list);
1612 case eenhENERGY_NSTEPS:
1613 do_cpt_step_err(xd, eenh_names[i], &enerhist->nsteps, list);
1615 case eenhENERGY_NSTEPS_SIM:
1616 do_cpt_step_err(xd, eenh_names[i], &enerhist->nsteps_sim, list);
1618 case eenhENERGY_DELTA_H_NN:
1621 if (!bRead && deltaH != nullptr)
1623 numDeltaH = deltaH->dh.size();
1625 do_cpt_int_err(xd, eenh_names[i], &numDeltaH, list);
1628 if (deltaH == nullptr)
1630 enerhist->deltaHForeignLambdas = std::make_unique<delta_h_history_t>();
1631 deltaH = enerhist->deltaHForeignLambdas.get();
1633 deltaH->dh.resize(numDeltaH);
1634 deltaH->start_lambda_set = FALSE;
1638 case eenhENERGY_DELTA_H_LIST:
1639 for (auto dh : deltaH->dh)
1641 ret = doVector<real>(xd, part, i, fflags, &dh, list);
1644 case eenhENERGY_DELTA_H_STARTTIME:
1645 ret = do_cpte_double(xd, part, i, fflags, &(deltaH->start_time), list);
1647 case eenhENERGY_DELTA_H_STARTLAMBDA:
1648 ret = do_cpte_double(xd, part, i, fflags, &(deltaH->start_lambda), list);
1652 "Unknown energy history entry %d\n"
1653 "You are probably reading a new checkpoint file with old code",
1659 if ((fflags & (1 << eenhENERGY_SUM)) && !(fflags & (1 << eenhENERGY_SUM_SIM)))
1661 /* Assume we have an old file format and copy sum to sum_sim */
1662 enerhist->ener_sum_sim = enerhist->ener_sum;
1665 if ((fflags & (1 << eenhENERGY_NSUM)) && !(fflags & (1 << eenhENERGY_NSTEPS)))
1667 /* Assume we have an old file format and copy nsum to nsteps */
1668 enerhist->nsteps = enerhist->nsum;
1670 if ((fflags & (1 << eenhENERGY_NSUM_SIM)) && !(fflags & (1 << eenhENERGY_NSTEPS_SIM)))
1672 /* Assume we have an old file format and copy nsum to nsteps */
1673 enerhist->nsteps_sim = enerhist->nsum_sim;
1679 static int doCptPullCoordHist(XDR* xd, PullCoordinateHistory* pullCoordHist, const StatePart part, FILE* list)
1684 flags |= ((1 << epullcoordh_VALUE_REF_SUM) | (1 << epullcoordh_VALUE_SUM)
1685 | (1 << epullcoordh_DR01_SUM) | (1 << epullcoordh_DR23_SUM)
1686 | (1 << epullcoordh_DR45_SUM) | (1 << epullcoordh_FSCAL_SUM));
1688 for (int i = 0; i < epullcoordh_NR && ret == 0; i++)
1692 case epullcoordh_VALUE_REF_SUM:
1693 ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->valueRef), list);
1695 case epullcoordh_VALUE_SUM:
1696 ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->value), list);
1698 case epullcoordh_DR01_SUM:
1699 for (int j = 0; j < DIM && ret == 0; j++)
1701 ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dr01[j]), list);
1704 case epullcoordh_DR23_SUM:
1705 for (int j = 0; j < DIM && ret == 0; j++)
1707 ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dr23[j]), list);
1710 case epullcoordh_DR45_SUM:
1711 for (int j = 0; j < DIM && ret == 0; j++)
1713 ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dr45[j]), list);
1716 case epullcoordh_FSCAL_SUM:
1717 ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->scalarForce), list);
1719 case epullcoordh_DYNAX_SUM:
1720 for (int j = 0; j < DIM && ret == 0; j++)
1722 ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dynaX[j]), list);
1731 static int doCptPullGroupHist(XDR* xd, PullGroupHistory* pullGroupHist, const StatePart part, FILE* list)
1736 flags |= ((1 << epullgrouph_X_SUM));
1738 for (int i = 0; i < epullgrouph_NR; i++)
1742 case epullgrouph_X_SUM:
1743 for (int j = 0; j < DIM && ret == 0; j++)
1745 ret = do_cpte_double(xd, part, i, flags, &(pullGroupHist->x[j]), list);
1755 static int doCptPullHist(XDR* xd, gmx_bool bRead, int fflags, PullHistory* pullHist, const StatePart part, FILE* list)
1758 int pullHistoryNumCoordinates = 0;
1759 int pullHistoryNumGroups = 0;
1761 /* Retain the number of terms in the sum and the number of coordinates (used for writing
1762 * average pull forces and coordinates) in the pull history, in temporary variables,
1763 * in case they cannot be read from the checkpoint, in order to have backward compatibility */
1766 pullHist->numValuesInXSum = 0;
1767 pullHist->numValuesInFSum = 0;
1769 else if (pullHist != nullptr)
1771 pullHistoryNumCoordinates = pullHist->pullCoordinateSums.size();
1772 pullHistoryNumGroups = pullHist->pullGroupSums.size();
1776 GMX_RELEASE_ASSERT(fflags == 0, "Without pull history, all flags should be off");
1779 for (int i = 0; (i < epullhNR && ret == 0); i++)
1781 if (fflags & (1 << i))
1785 case epullhPULL_NUMCOORDINATES:
1786 ret = do_cpte_int(xd, part, i, fflags, &pullHistoryNumCoordinates, list);
1788 case epullhPULL_NUMGROUPS:
1789 do_cpt_int_err(xd, eenh_names[i], &pullHistoryNumGroups, list);
1791 case epullhPULL_NUMVALUESINXSUM:
1792 do_cpt_int_err(xd, eenh_names[i], &pullHist->numValuesInXSum, list);
1794 case epullhPULL_NUMVALUESINFSUM:
1795 do_cpt_int_err(xd, eenh_names[i], &pullHist->numValuesInFSum, list);
1799 "Unknown pull history entry %d\n"
1800 "You are probably reading a new checkpoint file with old code",
1807 pullHist->pullCoordinateSums.resize(pullHistoryNumCoordinates);
1808 pullHist->pullGroupSums.resize(pullHistoryNumGroups);
1810 if (pullHist->numValuesInXSum > 0 || pullHist->numValuesInFSum > 0)
1812 for (size_t i = 0; i < pullHist->pullCoordinateSums.size() && ret == 0; i++)
1814 ret = doCptPullCoordHist(xd, &(pullHist->pullCoordinateSums[i]), part, list);
1816 for (size_t i = 0; i < pullHist->pullGroupSums.size() && ret == 0; i++)
1818 ret = doCptPullGroupHist(xd, &(pullHist->pullGroupSums[i]), part, list);
1825 static int do_cpt_df_hist(XDR* xd, int fflags, int nlambda, df_history_t** dfhistPtr, FILE* list)
1834 if (*dfhistPtr == nullptr)
1836 snew(*dfhistPtr, 1);
1837 (*dfhistPtr)->nlambda = nlambda;
1838 init_df_history(*dfhistPtr, nlambda);
1840 df_history_t* dfhist = *dfhistPtr;
1842 const StatePart part = StatePart::freeEnergyHistory;
1843 for (int i = 0; (i < edfhNR && ret == 0); i++)
1845 if (fflags & (1 << i))
1850 ret = do_cpte_bool(xd, part, i, fflags, &dfhist->bEquil, list);
1853 ret = do_cpte_ints(xd, part, i, fflags, nlambda, &dfhist->n_at_lam, list);
1856 ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->wl_histo, list);
1859 ret = do_cpte_real(xd, part, i, fflags, &dfhist->wl_delta, list);
1861 case edfhSUMWEIGHTS:
1862 ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->sum_weights, list);
1865 ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->sum_dg, list);
1868 ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->sum_minvar, list);
1871 ret = do_cpte_reals(xd, part, i, fflags, nlambda, &dfhist->sum_variance, list);
1874 ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->accum_p, list);
1877 ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->accum_m, list);
1880 ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->accum_p2, list);
1883 ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->accum_m2, list);
1886 ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->Tij, list);
1889 ret = do_cpte_nmatrix(xd, part, i, fflags, nlambda, dfhist->Tij_empirical, list);
1894 "Unknown df history entry %d\n"
1895 "You are probably reading a new checkpoint file with old code",
1905 /* This function stores the last whole configuration of the reference and
1906 * average structure in the .cpt file
1908 static int do_cpt_EDstate(XDR* xd, gmx_bool bRead, int nED, edsamhistory_t* EDstate, FILE* list)
1915 EDstate->bFromCpt = bRead;
1918 /* When reading, init_edsam has not been called yet,
1919 * so we have to allocate memory first. */
1922 snew(EDstate->nref, EDstate->nED);
1923 snew(EDstate->old_sref, EDstate->nED);
1924 snew(EDstate->nav, EDstate->nED);
1925 snew(EDstate->old_sav, EDstate->nED);
1928 /* Read/write the last whole conformation of SREF and SAV for each ED dataset (usually only one) */
1929 for (int i = 0; i < EDstate->nED; i++)
1933 /* Reference structure SREF */
1934 sprintf(buf, "ED%d # of atoms in reference structure", i + 1);
1935 do_cpt_int_err(xd, buf, &EDstate->nref[i], list);
1936 sprintf(buf, "ED%d x_ref", i + 1);
1939 snew(EDstate->old_sref[i], EDstate->nref[i]);
1940 do_cpt_n_rvecs_err(xd, buf, EDstate->nref[i], EDstate->old_sref[i], list);
1944 do_cpt_n_rvecs_err(xd, buf, EDstate->nref[i], EDstate->old_sref_p[i], list);
1947 /* Average structure SAV */
1948 sprintf(buf, "ED%d # of atoms in average structure", i + 1);
1949 do_cpt_int_err(xd, buf, &EDstate->nav[i], list);
1950 sprintf(buf, "ED%d x_av", i + 1);
1953 snew(EDstate->old_sav[i], EDstate->nav[i]);
1954 do_cpt_n_rvecs_err(xd, buf, EDstate->nav[i], EDstate->old_sav[i], list);
1958 do_cpt_n_rvecs_err(xd, buf, EDstate->nav[i], EDstate->old_sav_p[i], list);
1965 static int do_cpt_correlation_grid(XDR* xd,
1967 gmx_unused int fflags,
1968 gmx::CorrelationGridHistory* corrGrid,
1974 do_cpt_int_err(xd, eawhh_names[eawhh], &(corrGrid->numCorrelationTensors), list);
1975 do_cpt_int_err(xd, eawhh_names[eawhh], &(corrGrid->tensorSize), list);
1976 do_cpt_int_err(xd, eawhh_names[eawhh], &(corrGrid->blockDataListSize), list);
1980 initCorrelationGridHistory(corrGrid, corrGrid->numCorrelationTensors, corrGrid->tensorSize,
1981 corrGrid->blockDataListSize);
1984 for (gmx::CorrelationBlockDataHistory& blockData : corrGrid->blockDataBuffer)
1986 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumWeight), list);
1987 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumSquareWeight), list);
1988 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumWeightX), list);
1989 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockSumWeightY), list);
1990 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksSquareBlockWeight), list);
1991 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksBlockSquareWeight), list);
1992 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksBlockWeightBlockWeightX), list);
1993 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.sumOverBlocksBlockWeightBlockWeightY), list);
1994 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.blockLength), list);
1995 do_cpt_int_err(xd, eawhh_names[eawhh], &(blockData.previousBlockIndex), list);
1996 do_cpt_double_err(xd, eawhh_names[eawhh], &(blockData.correlationIntegral), list);
2002 static int do_cpt_awh_bias(XDR* xd, gmx_bool bRead, int fflags, gmx::AwhBiasHistory* biasHistory, FILE* list)
2006 gmx::AwhBiasStateHistory* state = &biasHistory->state;
2007 for (int i = 0; (i < eawhhNR && ret == 0); i++)
2009 if (fflags & (1 << i))
2013 case eawhhIN_INITIAL:
2014 do_cpt_bool_err(xd, eawhh_names[i], &state->in_initial, list);
2016 case eawhhEQUILIBRATEHISTOGRAM:
2017 do_cpt_bool_err(xd, eawhh_names[i], &state->equilibrateHistogram, list);
2020 do_cpt_double_err(xd, eawhh_names[i], &state->histSize, list);
2027 numPoints = biasHistory->pointState.size();
2029 do_cpt_int_err(xd, eawhh_names[i], &numPoints, list);
2032 biasHistory->pointState.resize(numPoints);
2036 case eawhhCOORDPOINT:
2037 for (auto& psh : biasHistory->pointState)
2039 do_cpt_double_err(xd, eawhh_names[i], &psh.target, list);
2040 do_cpt_double_err(xd, eawhh_names[i], &psh.free_energy, list);
2041 do_cpt_double_err(xd, eawhh_names[i], &psh.bias, list);
2042 do_cpt_double_err(xd, eawhh_names[i], &psh.weightsum_iteration, list);
2043 do_cpt_double_err(xd, eawhh_names[i], &psh.weightsum_covering, list);
2044 do_cpt_double_err(xd, eawhh_names[i], &psh.weightsum_tot, list);
2045 do_cpt_double_err(xd, eawhh_names[i], &psh.weightsum_ref, list);
2046 do_cpt_step_err(xd, eawhh_names[i], &psh.last_update_index, list);
2047 do_cpt_double_err(xd, eawhh_names[i], &psh.log_pmfsum, list);
2048 do_cpt_double_err(xd, eawhh_names[i], &psh.visits_iteration, list);
2049 do_cpt_double_err(xd, eawhh_names[i], &psh.visits_tot, list);
2052 case eawhhUMBRELLAGRIDPOINT:
2053 do_cpt_int_err(xd, eawhh_names[i], &(state->umbrellaGridpoint), list);
2055 case eawhhUPDATELIST:
2056 do_cpt_int_err(xd, eawhh_names[i], &(state->origin_index_updatelist), list);
2057 do_cpt_int_err(xd, eawhh_names[i], &(state->end_index_updatelist), list);
2059 case eawhhLOGSCALEDSAMPLEWEIGHT:
2060 do_cpt_double_err(xd, eawhh_names[i], &(state->logScaledSampleWeight), list);
2061 do_cpt_double_err(xd, eawhh_names[i], &(state->maxLogScaledSampleWeight), list);
2063 case eawhhNUMUPDATES:
2064 do_cpt_step_err(xd, eawhh_names[i], &(state->numUpdates), list);
2066 case eawhhFORCECORRELATIONGRID:
2067 ret = do_cpt_correlation_grid(xd, bRead, fflags,
2068 &biasHistory->forceCorrelationGrid, list, i);
2070 default: gmx_fatal(FARGS, "Unknown awh history entry %d\n", i);
2078 static int do_cpt_awh(XDR* xd, gmx_bool bRead, int fflags, gmx::AwhHistory* awhHistory, FILE* list)
2084 std::shared_ptr<gmx::AwhHistory> awhHistoryLocal;
2086 if (awhHistory == nullptr)
2088 GMX_RELEASE_ASSERT(bRead,
2089 "do_cpt_awh should not be called for writing without an AwhHistory");
2091 awhHistoryLocal = std::make_shared<gmx::AwhHistory>();
2092 awhHistory = awhHistoryLocal.get();
2095 /* To be able to read and write the AWH data several parameters determining the layout of the AWH structs need to be known
2096 (nbias, npoints, etc.). The best thing (?) would be to have these passed to this function. When writing to a checkpoint
2097 these parameters are available in awh_history (after calling init_awh_history). When reading from a checkpoint though, there
2098 is no initialized awh_history (it is initialized and set in this function). The AWH parameters have not always been read
2099 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
2100 when writing a checkpoint, also storing parameters needed for future reading of the checkpoint.
2102 Another issue is that some variables that AWH checkpoints don't have a registered enum and string (e.g. nbias below).
2103 One difficulty is the multilevel structure of the data which would need to be represented somehow. */
2108 numBias = awhHistory->bias.size();
2110 do_cpt_int_err(xd, "awh_nbias", &numBias, list);
2114 awhHistory->bias.resize(numBias);
2116 for (auto& bias : awhHistory->bias)
2118 ret = do_cpt_awh_bias(xd, bRead, fflags, &bias, list);
2124 do_cpt_double_err(xd, "awh_potential_offset", &awhHistory->potentialOffset, list);
2129 static void do_cpt_mdmodules(int fileVersion,
2130 t_fileio* checkpointFileHandle,
2131 const gmx::MdModulesNotifier& mdModulesNotifier)
2133 if (fileVersion >= cptv_MdModules)
2135 gmx::FileIOXdrSerializer serializer(checkpointFileHandle);
2136 gmx::KeyValueTreeObject mdModuleCheckpointParameterTree =
2137 gmx::deserializeKeyValueTree(&serializer);
2138 gmx::MdModulesCheckpointReadingDataOnMaster mdModuleCheckpointReadingDataOnMaster = {
2139 mdModuleCheckpointParameterTree, fileVersion
2141 mdModulesNotifier.checkpointingNotifications_.notify(mdModuleCheckpointReadingDataOnMaster);
2145 static int do_cpt_files(XDR* xd, gmx_bool bRead, std::vector<gmx_file_position_t>* outputfiles, FILE* list, int file_version)
2148 gmx_off_t mask = 0xFFFFFFFFL;
2149 int offset_high, offset_low;
2150 std::array<char, CPTSTRLEN> buf;
2151 GMX_RELEASE_ASSERT(outputfiles, "Must have valid outputfiles");
2153 // Ensure that reading pre-allocates outputfiles, while writing
2154 // writes what is already there.
2155 int nfiles = outputfiles->size();
2156 if (do_cpt_int(xd, "number of output files", &nfiles, list) != 0)
2162 outputfiles->resize(nfiles);
2165 for (auto& outputfile : *outputfiles)
2167 /* 64-bit XDR numbers are not portable, so it is stored as separate high/low fractions */
2170 do_cpt_string_err(xd, "output filename", buf, list);
2171 std::copy(std::begin(buf), std::end(buf), std::begin(outputfile.filename));
2173 if (do_cpt_int(xd, "file_offset_high", &offset_high, list) != 0)
2177 if (do_cpt_int(xd, "file_offset_low", &offset_low, list) != 0)
2181 outputfile.offset = (static_cast<gmx_off_t>(offset_high) << 32)
2182 | (static_cast<gmx_off_t>(offset_low) & mask);
2186 do_cpt_string_err(xd, "output filename", outputfile.filename, list);
2188 offset = outputfile.offset;
2196 offset_low = static_cast<int>(offset & mask);
2197 offset_high = static_cast<int>((offset >> 32) & mask);
2199 if (do_cpt_int(xd, "file_offset_high", &offset_high, list) != 0)
2203 if (do_cpt_int(xd, "file_offset_low", &offset_low, list) != 0)
2208 if (file_version >= 8)
2210 if (do_cpt_int(xd, "file_checksum_size", &outputfile.checksumSize, list) != 0)
2214 if (do_cpt_u_chars(xd, "file_checksum", outputfile.checksum.size(),
2215 outputfile.checksum.data(), list)
2223 outputfile.checksumSize = -1;
2229 void write_checkpoint_data(t_fileio* fp,
2230 CheckpointHeaderContents headerContents,
2234 ObservablesHistory* observablesHistory,
2235 const gmx::MdModulesNotifier& mdModulesNotifier,
2236 std::vector<gmx_file_position_t>* outputfiles,
2237 gmx::WriteCheckpointDataHolder* modularSimulatorCheckpointData)
2239 headerContents.flags_eks = 0;
2240 if (state->ekinstate.bUpToDate)
2242 headerContents.flags_eks = ((1 << eeksEKIN_N) | (1 << eeksEKINH) | (1 << eeksEKINF)
2243 | (1 << eeksEKINO) | (1 << eeksEKINSCALEF) | (1 << eeksEKINSCALEH)
2244 | (1 << eeksVSCALE) | (1 << eeksDEKINDL) | (1 << eeksMVCOS));
2246 headerContents.isModularSimulatorCheckpoint = !modularSimulatorCheckpointData->empty();
2248 energyhistory_t* enerhist = observablesHistory->energyHistory.get();
2249 headerContents.flags_enh = 0;
2250 if (enerhist != nullptr && (enerhist->nsum > 0 || enerhist->nsum_sim > 0))
2252 headerContents.flags_enh |=
2253 (1 << eenhENERGY_N) | (1 << eenhENERGY_NSTEPS) | (1 << eenhENERGY_NSTEPS_SIM);
2254 if (enerhist->nsum > 0)
2256 headerContents.flags_enh |=
2257 ((1 << eenhENERGY_AVER) | (1 << eenhENERGY_SUM) | (1 << eenhENERGY_NSUM));
2259 if (enerhist->nsum_sim > 0)
2261 headerContents.flags_enh |= ((1 << eenhENERGY_SUM_SIM) | (1 << eenhENERGY_NSUM_SIM));
2263 if (enerhist->deltaHForeignLambdas != nullptr)
2265 headerContents.flags_enh |=
2266 ((1 << eenhENERGY_DELTA_H_NN) | (1 << eenhENERGY_DELTA_H_LIST)
2267 | (1 << eenhENERGY_DELTA_H_STARTTIME) | (1 << eenhENERGY_DELTA_H_STARTLAMBDA));
2271 PullHistory* pullHist = observablesHistory->pullHistory.get();
2272 headerContents.flagsPullHistory = 0;
2273 if (pullHist != nullptr && (pullHist->numValuesInXSum > 0 || pullHist->numValuesInFSum > 0))
2275 headerContents.flagsPullHistory |= (1 << epullhPULL_NUMCOORDINATES);
2276 headerContents.flagsPullHistory |= ((1 << epullhPULL_NUMGROUPS) | (1 << epullhPULL_NUMVALUESINXSUM)
2277 | (1 << epullhPULL_NUMVALUESINFSUM));
2280 headerContents.flags_dfh = 0;
2283 headerContents.flags_dfh = ((1 << edfhBEQUIL) | (1 << edfhNATLAMBDA) | (1 << edfhSUMWEIGHTS)
2284 | (1 << edfhSUMDG) | (1 << edfhTIJ) | (1 << edfhTIJEMP));
2287 headerContents.flags_dfh |= ((1 << edfhWLDELTA) | (1 << edfhWLHISTO));
2289 if ((elamstats == elamstatsMINVAR) || (elamstats == elamstatsBARKER)
2290 || (elamstats == elamstatsMETROPOLIS))
2292 headerContents.flags_dfh |= ((1 << edfhACCUMP) | (1 << edfhACCUMM) | (1 << edfhACCUMP2)
2293 | (1 << edfhACCUMM2) | (1 << edfhSUMMINVAR) | (1 << edfhSUMVAR));
2297 headerContents.flags_awhh = 0;
2298 if (state->awhHistory != nullptr && !state->awhHistory->bias.empty())
2300 headerContents.flags_awhh |=
2301 ((1 << eawhhIN_INITIAL) | (1 << eawhhEQUILIBRATEHISTOGRAM) | (1 << eawhhHISTSIZE)
2302 | (1 << eawhhNPOINTS) | (1 << eawhhCOORDPOINT) | (1 << eawhhUMBRELLAGRIDPOINT)
2303 | (1 << eawhhUPDATELIST) | (1 << eawhhLOGSCALEDSAMPLEWEIGHT)
2304 | (1 << eawhhNUMUPDATES) | (1 << eawhhFORCECORRELATIONGRID));
2307 do_cpt_header(gmx_fio_getxdr(fp), FALSE, nullptr, &headerContents);
2309 if ((do_cpt_state(gmx_fio_getxdr(fp), state->flags, state, nullptr) < 0)
2310 || (do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents.flags_eks, &state->ekinstate, nullptr) < 0)
2311 || (do_cpt_enerhist(gmx_fio_getxdr(fp), FALSE, headerContents.flags_enh, enerhist, nullptr) < 0)
2312 || (doCptPullHist(gmx_fio_getxdr(fp), FALSE, headerContents.flagsPullHistory, pullHist,
2313 StatePart::pullHistory, nullptr)
2315 || (do_cpt_df_hist(gmx_fio_getxdr(fp), headerContents.flags_dfh, headerContents.nlambda,
2316 &state->dfhist, nullptr)
2318 || (do_cpt_EDstate(gmx_fio_getxdr(fp), FALSE, headerContents.nED,
2319 observablesHistory->edsamHistory.get(), nullptr)
2321 || (do_cpt_awh(gmx_fio_getxdr(fp), FALSE, headerContents.flags_awhh, state->awhHistory.get(), nullptr) < 0)
2322 || (do_cpt_swapstate(gmx_fio_getxdr(fp), FALSE, headerContents.eSwapCoords,
2323 observablesHistory->swapHistory.get(), nullptr)
2325 || (do_cpt_files(gmx_fio_getxdr(fp), FALSE, outputfiles, nullptr, headerContents.file_version) < 0))
2327 gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
2330 // Checkpointing MdModules
2332 gmx::KeyValueTreeBuilder builder;
2333 gmx::MdModulesWriteCheckpointData mdModulesWriteCheckpoint = { builder.rootObject(),
2334 headerContents.file_version };
2335 mdModulesNotifier.checkpointingNotifications_.notify(mdModulesWriteCheckpoint);
2336 auto tree = builder.build();
2337 gmx::FileIOXdrSerializer serializer(fp);
2338 gmx::serializeKeyValueTree(tree, &serializer);
2341 // Checkpointing modular simulator
2343 gmx::FileIOXdrSerializer serializer(fp);
2344 modularSimulatorCheckpointData->serialize(&serializer);
2347 do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
2349 /* Always FAH checkpoint immediately after a Gromacs checkpoint.
2351 * Note that it is critical that we save a FAH checkpoint directly
2352 * after writing a Gromacs checkpoint. If the program dies, either
2353 * by the machine powering off suddenly or the process being,
2354 * killed, FAH can recover files that have only appended data by
2355 * truncating them to the last recorded length. The Gromacs
2356 * checkpoint does not just append data, it is fully rewritten each
2357 * time so a crash between moving the new Gromacs checkpoint file in
2358 * to place and writing a FAH checkpoint is not recoverable. Thus
2359 * the time between these operations must be kept as short a
2366 static void check_int(FILE* fplog, const char* type, int p, int f, gmx_bool* mm)
2368 bool foundMismatch = (p != f);
2376 fprintf(fplog, " %s mismatch,\n", type);
2377 fprintf(fplog, " current program: %d\n", p);
2378 fprintf(fplog, " checkpoint file: %d\n", f);
2379 fprintf(fplog, "\n");
2383 static void check_string(FILE* fplog, const char* type, const char* p, const char* f, gmx_bool* mm)
2385 bool foundMismatch = (std::strcmp(p, f) != 0);
2393 fprintf(fplog, " %s mismatch,\n", type);
2394 fprintf(fplog, " current program: %s\n", p);
2395 fprintf(fplog, " checkpoint file: %s\n", f);
2396 fprintf(fplog, "\n");
2400 static void check_match(FILE* fplog,
2401 const t_commrec* cr,
2403 const CheckpointHeaderContents& headerContents,
2404 gmx_bool reproducibilityRequested)
2406 /* Note that this check_string on the version will also print a message
2407 * when only the minor version differs. But we only print a warning
2408 * message further down with reproducibilityRequested=TRUE.
2410 gmx_bool versionDiffers = FALSE;
2411 check_string(fplog, "Version", gmx_version(), headerContents.version, &versionDiffers);
2413 gmx_bool precisionDiffers = FALSE;
2414 check_int(fplog, "Double prec.", GMX_DOUBLE, headerContents.double_prec, &precisionDiffers);
2415 if (precisionDiffers)
2417 const char msg_precision_difference[] =
2418 "You are continuing a simulation with a different precision. Not matching\n"
2419 "single/double precision will lead to precision or performance loss.\n";
2422 fprintf(fplog, "%s\n", msg_precision_difference);
2426 gmx_bool mm = (versionDiffers || precisionDiffers);
2428 if (reproducibilityRequested)
2430 check_string(fplog, "Program name", gmx::getProgramContext().fullBinaryPath(),
2431 headerContents.fprog, &mm);
2433 check_int(fplog, "#ranks", cr->nnodes, headerContents.nnodes, &mm);
2436 if (cr->sizeOfDefaultCommunicator > 1 && reproducibilityRequested)
2438 // TODO: These checks are incorrect (see redmine #3309)
2439 check_int(fplog, "#PME-ranks", cr->npmenodes, headerContents.npme, &mm);
2441 int npp = cr->sizeOfDefaultCommunicator;
2442 if (cr->npmenodes >= 0)
2444 npp -= cr->npmenodes;
2446 int npp_f = headerContents.nnodes;
2447 if (headerContents.npme >= 0)
2449 npp_f -= headerContents.npme;
2453 check_int(fplog, "#DD-cells[x]", dd_nc[XX], headerContents.dd_nc[XX], &mm);
2454 check_int(fplog, "#DD-cells[y]", dd_nc[YY], headerContents.dd_nc[YY], &mm);
2455 check_int(fplog, "#DD-cells[z]", dd_nc[ZZ], headerContents.dd_nc[ZZ], &mm);
2461 /* Gromacs should be able to continue from checkpoints between
2462 * different patch level versions, but we do not guarantee
2463 * compatibility between different major/minor versions - check this.
2467 sscanf(gmx_version(), "%5d", &gmx_major);
2468 int ret = sscanf(headerContents.version, "%5d", &cpt_major);
2469 gmx_bool majorVersionDiffers = (ret < 1 || gmx_major != cpt_major);
2471 const char msg_major_version_difference[] =
2472 "The current GROMACS major version is not identical to the one that\n"
2473 "generated the checkpoint file. In principle GROMACS does not support\n"
2474 "continuation from checkpoints between different versions, so we advise\n"
2475 "against this. If you still want to try your luck we recommend that you use\n"
2476 "the -noappend flag to keep your output files from the two versions separate.\n"
2477 "This might also work around errors where the output fields in the energy\n"
2478 "file have changed between the different versions.\n";
2480 const char msg_mismatch_notice[] =
2481 "GROMACS patchlevel, binary or parallel settings differ from previous run.\n"
2482 "Continuation is exact, but not guaranteed to be binary identical.\n";
2484 if (majorVersionDiffers)
2488 fprintf(fplog, "%s\n", msg_major_version_difference);
2491 else if (reproducibilityRequested)
2493 /* Major & minor versions match at least, but something is different. */
2496 fprintf(fplog, "%s\n", msg_mismatch_notice);
2502 static void read_checkpoint(const char* fn,
2504 const t_commrec* cr,
2507 int* init_fep_state,
2508 CheckpointHeaderContents* headerContents,
2510 ObservablesHistory* observablesHistory,
2511 gmx_bool reproducibilityRequested,
2512 const gmx::MdModulesNotifier& mdModulesNotifier,
2513 gmx::ReadCheckpointDataHolder* modularSimulatorCheckpointData,
2514 bool useModularSimulator)
2517 char buf[STEPSTRSIZE];
2520 fp = gmx_fio_open(fn, "r");
2521 do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, headerContents);
2523 // If we are appending, then we don't want write to the open log
2524 // file because we still need to compute a checksum for it. In
2525 // that case, the filehandle will be nullptr. Otherwise, we report
2526 // to the new log file about the checkpoint file that we are
2528 FILE* fplog = gmx_fio_getfp(logfio);
2531 fprintf(fplog, "\n");
2532 fprintf(fplog, "Reading checkpoint file %s\n", fn);
2533 fprintf(fplog, " file generated by: %s\n", headerContents->fprog);
2534 fprintf(fplog, " file generated at: %s\n", headerContents->ftime);
2535 fprintf(fplog, " GROMACS double prec.: %d\n", headerContents->double_prec);
2536 fprintf(fplog, " simulation part #: %d\n", headerContents->simulation_part);
2537 fprintf(fplog, " step: %s\n", gmx_step_str(headerContents->step, buf));
2538 fprintf(fplog, " time: %f\n", headerContents->t);
2539 fprintf(fplog, "\n");
2542 if (headerContents->natoms != state->natoms)
2545 "Checkpoint file is for a system of %d atoms, while the current system consists "
2547 headerContents->natoms, state->natoms);
2549 if (headerContents->ngtc != state->ngtc)
2552 "Checkpoint file is for a system of %d T-coupling groups, while the current "
2553 "system consists of %d T-coupling groups",
2554 headerContents->ngtc, state->ngtc);
2556 if (headerContents->nnhpres != state->nnhpres)
2559 "Checkpoint file is for a system of %d NH-pressure-coupling variables, while the "
2560 "current system consists of %d NH-pressure-coupling variables",
2561 headerContents->nnhpres, state->nnhpres);
2564 int nlambdaHistory = (state->dfhist ? state->dfhist->nlambda : 0);
2565 if (headerContents->nlambda != nlambdaHistory)
2568 "Checkpoint file is for a system with %d lambda states, while the current system "
2569 "consists of %d lambda states",
2570 headerContents->nlambda, nlambdaHistory);
2573 init_gtc_state(state, state->ngtc, state->nnhpres,
2574 headerContents->nhchainlength); /* need to keep this here to keep the tpr format working */
2575 /* write over whatever was read; we use the number of Nose-Hoover chains from the checkpoint */
2577 if (headerContents->eIntegrator != eIntegrator)
2580 "Cannot change integrator during a checkpoint restart. Perhaps you should make a "
2581 "new .tpr with grompp -f new.mdp -t %s",
2585 // For modular simulator, no state object is populated, so we cannot do this check here!
2586 if (headerContents->flags_state != state->flags && !useModularSimulator)
2589 "Cannot change a simulation algorithm during a checkpoint restart. Perhaps you "
2590 "should make a new .tpr with grompp -f new.mdp -t %s",
2594 GMX_ASSERT(!(headerContents->isModularSimulatorCheckpoint && !useModularSimulator),
2595 "Checkpoint file was written by modular simulator, but the current simulation uses "
2596 "the legacy simulator.");
2597 GMX_ASSERT(!(!headerContents->isModularSimulatorCheckpoint && useModularSimulator),
2598 "Checkpoint file was written by legacy simulator, but the current simulation uses "
2599 "the modular simulator.");
2603 check_match(fplog, cr, dd_nc, *headerContents, reproducibilityRequested);
2606 ret = do_cpt_state(gmx_fio_getxdr(fp), headerContents->flags_state, state, nullptr);
2607 *init_fep_state = state->fep_state; /* there should be a better way to do this than setting it
2608 here. Investigate for 5.0. */
2613 ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents->flags_eks, &state->ekinstate, nullptr);
2618 state->ekinstate.hasReadEkinState = (((headerContents->flags_eks & (1 << eeksEKINH)) != 0)
2619 || ((headerContents->flags_eks & (1 << eeksEKINF)) != 0)
2620 || ((headerContents->flags_eks & (1 << eeksEKINO)) != 0)
2621 || (((headerContents->flags_eks & (1 << eeksEKINSCALEF))
2622 | (headerContents->flags_eks & (1 << eeksEKINSCALEH))
2623 | (headerContents->flags_eks & (1 << eeksVSCALE)))
2626 if (headerContents->flags_enh && observablesHistory->energyHistory == nullptr)
2628 observablesHistory->energyHistory = std::make_unique<energyhistory_t>();
2630 ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE, headerContents->flags_enh,
2631 observablesHistory->energyHistory.get(), nullptr);
2637 if (headerContents->flagsPullHistory)
2639 if (observablesHistory->pullHistory == nullptr)
2641 observablesHistory->pullHistory = std::make_unique<PullHistory>();
2643 ret = doCptPullHist(gmx_fio_getxdr(fp), TRUE, headerContents->flagsPullHistory,
2644 observablesHistory->pullHistory.get(), StatePart::pullHistory, nullptr);
2651 if (headerContents->file_version < 6)
2654 "Continuing from checkpoint files written before GROMACS 4.5 is not supported");
2657 ret = do_cpt_df_hist(gmx_fio_getxdr(fp), headerContents->flags_dfh, headerContents->nlambda,
2658 &state->dfhist, nullptr);
2664 if (headerContents->nED > 0 && observablesHistory->edsamHistory == nullptr)
2666 observablesHistory->edsamHistory = std::make_unique<edsamhistory_t>(edsamhistory_t{});
2668 ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, headerContents->nED,
2669 observablesHistory->edsamHistory.get(), nullptr);
2675 if (headerContents->flags_awhh != 0 && state->awhHistory == nullptr)
2677 state->awhHistory = std::make_shared<gmx::AwhHistory>();
2679 ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE, headerContents->flags_awhh, state->awhHistory.get(), nullptr);
2685 if (headerContents->eSwapCoords != eswapNO && observablesHistory->swapHistory == nullptr)
2687 observablesHistory->swapHistory = std::make_unique<swaphistory_t>(swaphistory_t{});
2689 ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents->eSwapCoords,
2690 observablesHistory->swapHistory.get(), nullptr);
2696 std::vector<gmx_file_position_t> outputfiles;
2697 ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, nullptr, headerContents->file_version);
2702 do_cpt_mdmodules(headerContents->file_version, fp, mdModulesNotifier);
2703 if (headerContents->file_version >= cptv_ModularSimulator)
2705 gmx::FileIOXdrSerializer serializer(fp);
2706 modularSimulatorCheckpointData->deserialize(&serializer);
2708 ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents->file_version);
2713 if (gmx_fio_close(fp) != 0)
2715 gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
2720 void load_checkpoint(const char* fn,
2722 const t_commrec* cr,
2726 ObservablesHistory* observablesHistory,
2727 gmx_bool reproducibilityRequested,
2728 const gmx::MdModulesNotifier& mdModulesNotifier,
2729 gmx::ReadCheckpointDataHolder* modularSimulatorCheckpointData,
2730 bool useModularSimulator)
2732 CheckpointHeaderContents headerContents;
2735 /* Read the state from the checkpoint file */
2736 read_checkpoint(fn, logfio, cr, dd_nc, ir->eI, &(ir->fepvals->init_fep_state),
2737 &headerContents, state, observablesHistory, reproducibilityRequested,
2738 mdModulesNotifier, modularSimulatorCheckpointData, useModularSimulator);
2742 gmx_bcast(sizeof(headerContents.step), &headerContents.step, cr->mpiDefaultCommunicator);
2743 gmx::MdModulesCheckpointReadingBroadcast broadcastCheckPointData = {
2744 cr->mpiDefaultCommunicator, PAR(cr), headerContents.file_version
2746 mdModulesNotifier.checkpointingNotifications_.notify(broadcastCheckPointData);
2748 ir->bContinuation = TRUE;
2749 if (ir->nsteps >= 0)
2751 // TODO Should the following condition be <=? Currently if you
2752 // pass a checkpoint written by an normal completion to a restart,
2753 // mdrun will read all input, does some work but no steps, and
2754 // write successful output. But perhaps that is not desirable.
2755 // Note that we do not intend to support the use of mdrun
2756 // -nsteps to circumvent this condition.
2757 if (ir->nsteps + ir->init_step < headerContents.step)
2759 char buf[STEPSTRSIZE];
2760 std::string message =
2761 gmx::formatString("The input requested %s steps, ", gmx_step_str(ir->nsteps, buf));
2762 if (ir->init_step > 0)
2764 message += gmx::formatString("starting from step %s, ", gmx_step_str(ir->init_step, buf));
2766 message += gmx::formatString(
2767 "however the checkpoint "
2768 "file has already reached step %s. The simulation will not "
2769 "proceed, because either your simulation is already complete, "
2770 "or your combination of input files don't match.",
2771 gmx_step_str(headerContents.step, buf));
2772 gmx_fatal(FARGS, "%s", message.c_str());
2774 ir->nsteps += ir->init_step - headerContents.step;
2776 ir->init_step = headerContents.step;
2777 ir->simulation_part = headerContents.simulation_part + 1;
2780 void read_checkpoint_part_and_step(const char* filename, int* simulation_part, int64_t* step)
2784 if (filename == nullptr || !gmx_fexist(filename) || ((fp = gmx_fio_open(filename, "r")) == nullptr))
2786 *simulation_part = 0;
2791 CheckpointHeaderContents headerContents;
2792 do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, &headerContents);
2794 *simulation_part = headerContents.simulation_part;
2795 *step = headerContents.step;
2798 static CheckpointHeaderContents read_checkpoint_data(t_fileio* fp,
2800 std::vector<gmx_file_position_t>* outputfiles)
2802 CheckpointHeaderContents headerContents;
2803 do_cpt_header(gmx_fio_getxdr(fp), TRUE, nullptr, &headerContents);
2804 state->natoms = headerContents.natoms;
2805 state->ngtc = headerContents.ngtc;
2806 state->nnhpres = headerContents.nnhpres;
2807 state->nhchainlength = headerContents.nhchainlength;
2808 state->flags = headerContents.flags_state;
2809 int ret = do_cpt_state(gmx_fio_getxdr(fp), state->flags, state, nullptr);
2814 ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents.flags_eks, &state->ekinstate, nullptr);
2820 energyhistory_t enerhist;
2821 ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE, headerContents.flags_enh, &enerhist, nullptr);
2826 PullHistory pullHist = {};
2827 ret = doCptPullHist(gmx_fio_getxdr(fp), TRUE, headerContents.flagsPullHistory, &pullHist,
2828 StatePart::pullHistory, nullptr);
2834 ret = do_cpt_df_hist(gmx_fio_getxdr(fp), headerContents.flags_dfh, headerContents.nlambda,
2835 &state->dfhist, nullptr);
2841 edsamhistory_t edsamhist = {};
2842 ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, headerContents.nED, &edsamhist, nullptr);
2848 ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE, headerContents.flags_awhh, state->awhHistory.get(), nullptr);
2854 swaphistory_t swaphist = {};
2855 ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents.eSwapCoords, &swaphist, nullptr);
2861 ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, outputfiles, nullptr, headerContents.file_version);
2867 gmx::MdModulesNotifier mdModuleNotifier;
2868 do_cpt_mdmodules(headerContents.file_version, fp, mdModuleNotifier);
2869 if (headerContents.file_version >= cptv_ModularSimulator)
2871 // In the scope of the current function, we can just throw away the content
2872 // of the modular checkpoint, but we need to read it to move the file pointer
2873 gmx::FileIOXdrSerializer serializer(fp);
2874 gmx::ReadCheckpointDataHolder modularSimulatorCheckpointData;
2875 modularSimulatorCheckpointData.deserialize(&serializer);
2877 ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
2882 return headerContents;
2885 void read_checkpoint_trxframe(t_fileio* fp, t_trxframe* fr)
2888 std::vector<gmx_file_position_t> outputfiles;
2889 CheckpointHeaderContents headerContents = read_checkpoint_data(fp, &state, &outputfiles);
2891 fr->natoms = state.natoms;
2893 fr->step = int64_to_int(headerContents.step, "conversion of checkpoint to trajectory");
2895 fr->time = headerContents.t;
2897 fr->lambda = state.lambda[efptFEP];
2898 fr->fep_state = state.fep_state;
2900 fr->bX = ((state.flags & (1 << estX)) != 0);
2903 fr->x = makeRvecArray(state.x, state.natoms);
2905 fr->bV = ((state.flags & (1 << estV)) != 0);
2908 fr->v = makeRvecArray(state.v, state.natoms);
2911 fr->bBox = ((state.flags & (1 << estBOX)) != 0);
2914 copy_mat(state.box, fr->box);
2918 void list_checkpoint(const char* fn, FILE* out)
2925 fp = gmx_fio_open(fn, "r");
2926 CheckpointHeaderContents headerContents;
2927 do_cpt_header(gmx_fio_getxdr(fp), TRUE, out, &headerContents);
2928 state.natoms = headerContents.natoms;
2929 state.ngtc = headerContents.ngtc;
2930 state.nnhpres = headerContents.nnhpres;
2931 state.nhchainlength = headerContents.nhchainlength;
2932 state.flags = headerContents.flags_state;
2933 ret = do_cpt_state(gmx_fio_getxdr(fp), state.flags, &state, out);
2938 ret = do_cpt_ekinstate(gmx_fio_getxdr(fp), headerContents.flags_eks, &state.ekinstate, out);
2944 energyhistory_t enerhist;
2945 ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE, headerContents.flags_enh, &enerhist, out);
2949 PullHistory pullHist = {};
2950 ret = doCptPullHist(gmx_fio_getxdr(fp), TRUE, headerContents.flagsPullHistory, &pullHist,
2951 StatePart::pullHistory, out);
2956 ret = do_cpt_df_hist(gmx_fio_getxdr(fp), headerContents.flags_dfh, headerContents.nlambda,
2957 &state.dfhist, out);
2962 edsamhistory_t edsamhist = {};
2963 ret = do_cpt_EDstate(gmx_fio_getxdr(fp), TRUE, headerContents.nED, &edsamhist, out);
2968 ret = do_cpt_awh(gmx_fio_getxdr(fp), TRUE, headerContents.flags_awhh, state.awhHistory.get(), out);
2973 swaphistory_t swaphist = {};
2974 ret = do_cpt_swapstate(gmx_fio_getxdr(fp), TRUE, headerContents.eSwapCoords, &swaphist, out);
2979 std::vector<gmx_file_position_t> outputfiles;
2980 ret = do_cpt_files(gmx_fio_getxdr(fp), TRUE, &outputfiles, out, headerContents.file_version);
2985 ret = do_cpt_footer(gmx_fio_getxdr(fp), headerContents.file_version);
2992 if (gmx_fio_close(fp) != 0)
2994 gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
2998 /* This routine cannot print tons of data, since it is called before the log file is opened. */
2999 CheckpointHeaderContents read_checkpoint_simulation_part_and_filenames(t_fileio* fp,
3000 std::vector<gmx_file_position_t>* outputfiles)
3003 CheckpointHeaderContents headerContents = read_checkpoint_data(fp, &state, outputfiles);
3004 if (gmx_fio_close(fp) != 0)
3006 gmx_file("Cannot read/write checkpoint; corrupt file, or maybe you are out of disk space?");
3008 return headerContents;