src/gromacs/fileio/tpxio.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 /* This file is completely threadsafe - keep it that way! */
  41
  42 #include "gromacs/fileio/tpxio.h"
  43
  44 #include <cstdio>
  45 #include <cstdlib>
  46 #include <cstring>
  47
  48 #include <algorithm>
  49 #include <memory>
  50 #include <vector>
  51
  52 #include "gromacs/applied_forces/awh/read_params.h"
  53 #include "gromacs/fileio/filetypes.h"
  54 #include "gromacs/fileio/gmxfio.h"
  55 #include "gromacs/fileio/gmxfio_xdr.h"
  56 #include "gromacs/math/units.h"
  57 #include "gromacs/math/vec.h"
  58 #include "gromacs/mdtypes/awh_history.h"
  59 #include "gromacs/mdtypes/awh_params.h"
  60 #include "gromacs/mdtypes/inputrec.h"
  61 #include "gromacs/mdtypes/md_enums.h"
  62 #include "gromacs/mdtypes/multipletimestepping.h"
  63 #include "gromacs/mdtypes/pull_params.h"
  64 #include "gromacs/mdtypes/state.h"
  65 #include "gromacs/pbcutil/boxutilities.h"
  66 #include "gromacs/pbcutil/pbc.h"
  67 #include "gromacs/topology/block.h"
  68 #include "gromacs/topology/ifunc.h"
  69 #include "gromacs/topology/mtop_util.h"
  70 #include "gromacs/topology/symtab.h"
  71 #include "gromacs/topology/topology.h"
  72 #include "gromacs/utility/arraysize.h"
  73 #include "gromacs/utility/baseversion.h"
  74 #include "gromacs/utility/cstringutil.h"
  75 #include "gromacs/utility/fatalerror.h"
  76 #include "gromacs/utility/futil.h"
  77 #include "gromacs/utility/gmxassert.h"
  78 #include "gromacs/utility/inmemoryserializer.h"
  79 #include "gromacs/utility/iserializer.h"
  80 #include "gromacs/utility/keyvaluetreebuilder.h"
  81 #include "gromacs/utility/keyvaluetreeserializer.h"
  82 #include "gromacs/utility/smalloc.h"
  83 #include "gromacs/utility/snprintf.h"
  84 #include "gromacs/utility/txtdump.h"
  85
  86 #define TPX_TAG_RELEASE "release"
  87
  88 /*! \brief Tag string for the file format written to run input files
  89  * written by this version of the code.
  90  *
  91  * Change this if you want to change the run input format in a feature
  92  * branch. This ensures that there will not be different run input
  93  * formats around which cannot be distinguished, while not causing
  94  * problems rebasing the feature branch onto upstream changes. When
  95  * merging with mainstream GROMACS, set this tag string back to
  96  * TPX_TAG_RELEASE, and instead add an element to tpxv.
  97  */
  98 static const std::string tpx_tag = TPX_TAG_RELEASE;
  99
 100 /*! \brief Enum of values that describe the contents of a tpr file
 101  * whose format matches a version number
 102  *
 103  * The enum helps the code be more self-documenting and ensure merges
 104  * do not silently resolve when two patches make the same bump. When
 105  * adding new functionality, add a new element just above tpxv_Count
 106  * in this enumeration, and write code below that does the right thing
 107  * according to the value of file_version.
 108  */
 109 enum tpxv
 110 {
 111     tpxv_ComputationalElectrophysiology =
 112             96, /**< support for ion/water position swaps (computational electrophysiology) */
 113     tpxv_Use64BitRandomSeed, /**< change ld_seed from int to int64_t */
 114     tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, /**< potentials for supporting coarse-grained force fields */
 115     tpxv_InteractiveMolecularDynamics, /**< interactive molecular dynamics (IMD) */
 116     tpxv_RemoveObsoleteParameters1,    /**< remove optimize_fft, dihre_fc, nstcheckpoint */
 117     tpxv_PullCoordTypeGeom,            /**< add pull type and geometry per group and flat-bottom */
 118     tpxv_PullGeomDirRel,               /**< add pull geometry direction-relative */
 119     tpxv_IntermolecularBondeds, /**< permit inter-molecular bonded interactions in the topology */
 120     tpxv_CompElWithSwapLayerOffset, /**< added parameters for improved CompEl setups */
 121     tpxv_CompElPolyatomicIonsAndMultipleIonTypes, /**< CompEl now can handle polyatomic ions and more than two types of ions */
 122     tpxv_RemoveAdress,                            /**< removed support for AdResS */
 123     tpxv_PullCoordNGroup,               /**< add ngroup to pull coord */
 124     tpxv_RemoveTwinRange,               /**< removed support for twin-range interactions */
 125     tpxv_ReplacePullPrintCOM12,         /**< Replaced print-com-1, 2 with pull-print-com */
 126     tpxv_PullExternalPotential,         /**< Added pull type external potential */
 127     tpxv_GenericParamsForElectricField, /**< Introduced KeyValueTree and moved electric field parameters */
 128     tpxv_AcceleratedWeightHistogram, /**< sampling with accelerated weight histogram method (AWH) */
 129     tpxv_RemoveImplicitSolvation,    /**< removed support for implicit solvation */
 130     tpxv_PullPrevStepCOMAsReference, /**< Enabled using the COM of the pull group of the last frame as reference for PBC */
 131     tpxv_MimicQMMM,   /**< Introduced support for MiMiC QM/MM interface */
 132     tpxv_PullAverage, /**< Added possibility to output average pull force and position */
 133     tpxv_GenericInternalParameters, /**< Added internal parameters for mdrun modules*/
 134     tpxv_VSite2FD,                  /**< Added 2FD type virtual site */
 135     tpxv_AddSizeField, /**< Added field with information about the size of the serialized tpr file in bytes, excluding the header */
 136     tpxv_StoreNonBondedInteractionExclusionGroup, /**< Store the non bonded interaction exclusion group in the topology */
 137     tpxv_VSite1,                                  /**< Added 1 type virtual site */
 138     tpxv_MTS,                                     /**< Added multiple time stepping */
 139     tpxv_RemovedConstantAcceleration, /**< Removed support for constant acceleration NEMD. */
 140     tpxv_TransformationPullCoord,     /**< Support for transformation pull coordinates */
 141     tpxv_SoftcoreGapsys,              /**< Added gapsys softcore function */
 142     tpxv_ReaddedConstantAcceleration, /**< Re-added support for constant acceleration NEMD. */
 143     tpxv_Count                        /**< the total number of tpxv versions */
 144 };
 145
 146 /*! \brief Version number of the file format written to run input
 147  * files by this version of the code.
 148  *
 149  * The tpx_version increases whenever the file format in the main
 150  * development branch changes, due to an extension of the tpxv enum above.
 151  * Backward compatibility for reading old run input files is maintained
 152  * by checking this version number against that of the file and then using
 153  * the correct code path.
 154  *
 155  * When developing a feature branch that needs to change the run input
 156  * file format, change tpx_tag instead. */
 157 static const int tpx_version = tpxv_Count - 1;
 158
 159
 160 /*! \brief
 161  * Enum keeping track of incompatible changes for older TPR versions.
 162  *
 163  * The enum should be updated with a new field when editing the TOPOLOGY
 164  * or HEADER of the tpx format. In particular, updating ftupd or
 165  * changing the fields of TprHeaderVersion often trigger such needs.
 166  *
 167  * This way we can maintain forward compatibility too for all analysis tools
 168  * and/or external programs that only need to know the atom/residue names,
 169  * charges, and bond connectivity.
 170  *
 171  * It first appeared in tpx version 26, when I also moved the inputrecord
 172  * to the end of the tpx file, so we can just skip it if we only
 173  * want the topology.
 174  *
 175  * In particular, it must be increased when adding new elements to
 176  * ftupd, so that old code can read new .tpr files.
 177  */
 178 enum class TpxGeneration : int
 179 {
 180     Initial = 26, //! First version is 26
 181     AddSizeField, //! TPR header modified for writing as a block.
 182     AddVSite1,    //! ftupd changed to include VSite1 type.
 183     Count         //! Number of entries.
 184 };
 185
 186 //! Value of Current TPR generation.
 187 static const int tpx_generation = static_cast<int>(TpxGeneration::Count) - 1;
 188
 189 /* This number should be the most recent backwards incompatible version
 190  * I.e., if this number is 9, we cannot read tpx version 9 with this code.
 191  */
 192 static const int tpx_incompatible_version = 57; // GMX4.0 has version 58
 193
 194
 195 /* Struct used to maintain tpx compatibility when function types are added */
 196 typedef struct
 197 {
 198     int fvnr;  /* file version number in which the function type first appeared */
 199     int ftype; /* function type */
 200 } t_ftupd;
 201
 202 /*
 203  * TODO The following three lines make little sense, please clarify if
 204  * you've had to work out how ftupd works.
 205  *
 206  * The entries should be ordered in:
 207  * 1. ascending function type number
 208  * 2. ascending file version number
 209  *
 210  * Because we support reading of old .tpr file versions (even when
 211  * mdrun can no longer run the simulation), we need to be able to read
 212  * obsolete t_interaction_function types. Any data read from such
 213  * fields is discarded. Their names have _NOLONGERUSED appended to
 214  * them to make things clear.
 215  *
 216  * When adding to or making breaking changes to reading this struct,
 217  * update TpxGeneration.
 218  */
 219 static const t_ftupd ftupd[] = {
 220     { 70, F_RESTRBONDS },
 221     { tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, F_RESTRANGLES },
 222     { 76, F_LINEAR_ANGLES },
 223     { tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, F_RESTRDIHS },
 224     { tpxv_RestrictedBendingAndCombinedAngleTorsionPotentials, F_CBTDIHS },
 225     { 65, F_CMAP },
 226     { 60, F_GB12_NOLONGERUSED },
 227     { 61, F_GB13_NOLONGERUSED },
 228     { 61, F_GB14_NOLONGERUSED },
 229     { 72, F_GBPOL_NOLONGERUSED },
 230     { 72, F_NPSOLVATION_NOLONGERUSED },
 231     { 93, F_LJ_RECIP },
 232     { 76, F_ANHARM_POL },
 233     { 90, F_FBPOSRES },
 234     { tpxv_VSite1, F_VSITE1 },
 235     { tpxv_VSite2FD, F_VSITE2FD },
 236     { tpxv_GenericInternalParameters, F_DENSITYFITTING },
 237     { 69, F_VTEMP_NOLONGERUSED },
 238     { 66, F_PDISPCORR },
 239     { 79, F_DVDL_COUL },
 240     {
 241             79,
 242             F_DVDL_VDW,
 243     },
 244     {
 245             79,
 246             F_DVDL_BONDED,
 247     },
 248     { 79, F_DVDL_RESTRAINT },
 249     { 79, F_DVDL_TEMPERATURE },
 250 };
 251 #define NFTUPD asize(ftupd)
 252
 253 /**************************************************************
 254  *
 255  * Now the higer level routines that do io of the structures and arrays
 256  *
 257  **************************************************************/
 258 static void do_pullgrp_tpx_pre95(gmx::ISerializer* serializer, t_pull_group* pgrp, t_pull_coord* pcrd)
 259 {
 260     rvec tmp;
 261
 262     int numAtoms = pgrp->ind.size();
 263     serializer->doInt(&numAtoms);
 264     pgrp->ind.resize(numAtoms);
 265     serializer->doIntArray(pgrp->ind.data(), numAtoms);
 266     int numWeights = pgrp->weight.size();
 267     serializer->doInt(&numWeights);
 268     pgrp->weight.resize(numWeights);
 269     serializer->doRealArray(pgrp->weight.data(), numWeights);
 270     serializer->doInt(&pgrp->pbcatom);
 271     serializer->doRvec(&pcrd->vec.as_vec());
 272     clear_rvec(pcrd->origin);
 273     serializer->doRvec(&tmp);
 274     pcrd->init = tmp[0];
 275     serializer->doReal(&pcrd->rate);
 276     serializer->doReal(&pcrd->k);
 277     serializer->doReal(&pcrd->kB);
 278 }
 279
 280 static void do_pull_group(gmx::ISerializer* serializer, t_pull_group* pgrp)
 281 {
 282     int numAtoms = pgrp->ind.size();
 283     serializer->doInt(&numAtoms);
 284     pgrp->ind.resize(numAtoms);
 285     serializer->doIntArray(pgrp->ind.data(), numAtoms);
 286     int numWeights = pgrp->weight.size();
 287     serializer->doInt(&numWeights);
 288     pgrp->weight.resize(numWeights);
 289     serializer->doRealArray(pgrp->weight.data(), numWeights);
 290     serializer->doInt(&pgrp->pbcatom);
 291 }
 292
 293 static void do_pull_coord(gmx::ISerializer* serializer,
 294                           t_pull_coord*     pcrd,
 295                           int               file_version,
 296                           PullingAlgorithm  ePullOld,
 297                           PullGroupGeometry eGeomOld,
 298                           ivec              dimOld)
 299 {
 300     if (file_version >= tpxv_PullCoordNGroup)
 301     {
 302         serializer->doEnumAsInt(&pcrd->eType);
 303         if (file_version >= tpxv_PullExternalPotential)
 304         {
 305             if (pcrd->eType == PullingAlgorithm::External)
 306             {
 307                 serializer->doString(&pcrd->externalPotentialProvider);
 308             }
 309             else
 310             {
 311                 pcrd->externalPotentialProvider.clear();
 312             }
 313         }
 314         else
 315         {
 316             if (serializer->reading())
 317             {
 318                 pcrd->externalPotentialProvider.clear();
 319             }
 320         }
 321         /* Note that we try to support adding new geometries without
 322          * changing the tpx version. This requires checks when printing the
 323          * geometry string and a check and fatal_error in init_pull.
 324          */
 325         serializer->doEnumAsInt(&pcrd->eGeom);
 326         serializer->doInt(&pcrd->ngroup);
 327         if (pcrd->ngroup <= c_pullCoordNgroupMax)
 328         {
 329             serializer->doIntArray(pcrd->group.data(), pcrd->ngroup);
 330         }
 331         else
 332         {
 333             /* More groups in file than supported, this must be a new geometry
 334              * that is not supported by our current code. Since we will not
 335              * use the groups for this coord (checks in the pull and WHAM code
 336              * ensure this), we can ignore the groups and set ngroup=0.
 337              */
 338             int* dum;
 339             snew(dum, pcrd->ngroup);
 340             serializer->doIntArray(dum, pcrd->ngroup);
 341             sfree(dum);
 342
 343             pcrd->ngroup = 0;
 344         }
 345         serializer->doIvec(&pcrd->dim.as_vec());
 346         if (file_version >= tpxv_TransformationPullCoord)
 347         {
 348             serializer->doString(&pcrd->expression);
 349         }
 350         else
 351         {
 352             if (serializer->reading())
 353             {
 354                 pcrd->expression.clear();
 355             }
 356         }
 357     }
 358     else
 359     {
 360         pcrd->ngroup = 2;
 361         serializer->doInt(&pcrd->group[0]);
 362         serializer->doInt(&pcrd->group[1]);
 363         if (file_version >= tpxv_PullCoordTypeGeom)
 364         {
 365             pcrd->ngroup = (pcrd->eGeom == PullGroupGeometry::DirectionRelative ? 4 : 2);
 366             serializer->doEnumAsInt(&pcrd->eType);
 367             serializer->doEnumAsInt(&pcrd->eGeom);
 368             if (pcrd->ngroup == 4)
 369             {
 370                 serializer->doInt(&pcrd->group[2]);
 371                 serializer->doInt(&pcrd->group[3]);
 372             }
 373             serializer->doIvec(&pcrd->dim.as_vec());
 374         }
 375         else
 376         {
 377             pcrd->eType = ePullOld;
 378             pcrd->eGeom = eGeomOld;
 379             copy_ivec(dimOld, pcrd->dim);
 380         }
 381     }
 382     serializer->doRvec(&pcrd->origin.as_vec());
 383     serializer->doRvec(&pcrd->vec.as_vec());
 384     if (file_version >= tpxv_PullCoordTypeGeom)
 385     {
 386         serializer->doBool(&pcrd->bStart);
 387     }
 388     else
 389     {
 390         /* This parameter is only printed, but not actually used by mdrun */
 391         pcrd->bStart = FALSE;
 392     }
 393     serializer->doReal(&pcrd->init);
 394     serializer->doReal(&pcrd->rate);
 395     serializer->doReal(&pcrd->k);
 396     serializer->doReal(&pcrd->kB);
 397 }
 398
 399 static void do_expandedvals(gmx::ISerializer* serializer, t_expanded* expand, t_lambda* fepvals, int file_version)
 400 {
 401     int n_lambda = fepvals->n_lambda;
 402
 403     /* reset the lambda calculation window */
 404     fepvals->lambda_start_n = 0;
 405     fepvals->lambda_stop_n  = n_lambda;
 406     if (file_version >= 79)
 407     {
 408         if (n_lambda > 0)
 409         {
 410             expand->init_lambda_weights.resize(n_lambda);
 411             serializer->doRealArray(expand->init_lambda_weights.data(), n_lambda);
 412             serializer->doBool(&expand->bInit_weights);
 413         }
 414
 415         serializer->doInt(&expand->nstexpanded);
 416         serializer->doEnumAsInt(&expand->elmcmove);
 417         serializer->doEnumAsInt(&expand->elamstats);
 418         serializer->doInt(&expand->lmc_repeats);
 419         serializer->doInt(&expand->gibbsdeltalam);
 420         serializer->doInt(&expand->lmc_forced_nstart);
 421         serializer->doInt(&expand->lmc_seed);
 422         serializer->doReal(&expand->mc_temp);
 423         serializer->doBool(&expand->bSymmetrizedTMatrix);
 424         serializer->doInt(&expand->nstTij);
 425         serializer->doInt(&expand->minvarmin);
 426         serializer->doInt(&expand->c_range);
 427         serializer->doReal(&expand->wl_scale);
 428         serializer->doReal(&expand->wl_ratio);
 429         serializer->doReal(&expand->init_wl_delta);
 430         serializer->doBool(&expand->bWLoneovert);
 431         serializer->doEnumAsInt(&expand->elmceq);
 432         serializer->doInt(&expand->equil_steps);
 433         serializer->doInt(&expand->equil_samples);
 434         serializer->doInt(&expand->equil_n_at_lam);
 435         serializer->doReal(&expand->equil_wl_delta);
 436         serializer->doReal(&expand->equil_ratio);
 437     }
 438 }
 439
 440 static void do_simtempvals(gmx::ISerializer* serializer, t_simtemp* simtemp, int n_lambda, int file_version)
 441 {
 442     if (file_version >= 79)
 443     {
 444         serializer->doEnumAsInt(&simtemp->eSimTempScale);
 445         serializer->doReal(&simtemp->simtemp_high);
 446         serializer->doReal(&simtemp->simtemp_low);
 447         if (n_lambda > 0)
 448         {
 449             if (serializer->reading())
 450             {
 451                 simtemp->temperatures.resize(n_lambda);
 452             }
 453             serializer->doRealArray(simtemp->temperatures.data(), n_lambda);
 454         }
 455     }
 456 }
 457
 458 static void do_imd(gmx::ISerializer* serializer, t_IMD* imd)
 459 {
 460     serializer->doInt(&imd->nat);
 461     if (serializer->reading())
 462     {
 463         snew(imd->ind, imd->nat);
 464     }
 465     serializer->doIntArray(imd->ind, imd->nat);
 466 }
 467
 468 static void do_fepvals(gmx::ISerializer* serializer, t_lambda* fepvals, int file_version)
 469 {
 470     /* i is defined in the ndo_double macro; use g to iterate. */
 471     real rdum;
 472
 473     /* free energy values */
 474
 475     if (file_version >= 79)
 476     {
 477         serializer->doInt(&fepvals->init_fep_state);
 478         serializer->doDouble(&fepvals->init_lambda);
 479         serializer->doDouble(&fepvals->delta_lambda);
 480     }
 481     else if (file_version >= 59)
 482     {
 483         serializer->doDouble(&fepvals->init_lambda);
 484         serializer->doDouble(&fepvals->delta_lambda);
 485     }
 486     else
 487     {
 488         serializer->doReal(&rdum);
 489         fepvals->init_lambda = rdum;
 490         serializer->doReal(&rdum);
 491         fepvals->delta_lambda = rdum;
 492     }
 493     if (file_version >= 79)
 494     {
 495         serializer->doInt(&fepvals->n_lambda);
 496         for (auto g : keysOf(fepvals->all_lambda))
 497         {
 498             if (fepvals->n_lambda > 0)
 499             {
 500                 fepvals->all_lambda[g].resize(fepvals->n_lambda);
 501                 serializer->doDoubleArray(fepvals->all_lambda[g].data(), fepvals->n_lambda);
 502                 serializer->doBoolArray(fepvals->separate_dvdl.begin(), fepvals->separate_dvdl.size());
 503             }
 504             else if (fepvals->init_lambda >= 0)
 505             {
 506                 fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Fep] = TRUE;
 507             }
 508         }
 509     }
 510     else if (file_version >= 64)
 511     {
 512         serializer->doInt(&fepvals->n_lambda);
 513         if (serializer->reading())
 514         {
 515             /* still allocate the all_lambda array's contents. */
 516             for (auto g : keysOf(fepvals->all_lambda))
 517             {
 518                 if (fepvals->n_lambda > 0)
 519                 {
 520                     fepvals->all_lambda[g].resize(fepvals->n_lambda);
 521                 }
 522             }
 523         }
 524         serializer->doDoubleArray(fepvals->all_lambda[FreeEnergyPerturbationCouplingType::Fep].data(),
 525                                   fepvals->n_lambda);
 526         if (fepvals->init_lambda >= 0)
 527         {
 528             fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Fep] = TRUE;
 529
 530             if (serializer->reading())
 531             {
 532                 /* copy the contents of the efptFEP lambda component to all
 533                    the other components */
 534                 for (auto g : keysOf(fepvals->all_lambda))
 535                 {
 536                     for (int h = 0; h < fepvals->n_lambda; h++)
 537                     {
 538                         if (g != FreeEnergyPerturbationCouplingType::Fep)
 539                         {
 540                             fepvals->all_lambda[g][h] =
 541                                     fepvals->all_lambda[FreeEnergyPerturbationCouplingType::Fep][h];
 542                         }
 543                     }
 544                 }
 545             }
 546         }
 547     }
 548     else
 549     {
 550         fepvals->n_lambda = 0;
 551         if (fepvals->init_lambda >= 0)
 552         {
 553             fepvals->separate_dvdl[FreeEnergyPerturbationCouplingType::Fep] = TRUE;
 554         }
 555     }
 556     serializer->doReal(&fepvals->sc_alpha);
 557     serializer->doInt(&fepvals->sc_power);
 558     if (file_version >= 79)
 559     {
 560         serializer->doReal(&fepvals->sc_r_power);
 561     }
 562     else
 563     {
 564         fepvals->sc_r_power = 6.0;
 565     }
 566     if (fepvals->sc_r_power != 6.0)
 567     {
 568         gmx_fatal(FARGS, "sc-r-power=48 is no longer supported");
 569     }
 570     serializer->doReal(&fepvals->sc_sigma);
 571     if (serializer->reading())
 572     {
 573         if (file_version >= 71)
 574         {
 575             fepvals->sc_sigma_min = fepvals->sc_sigma;
 576         }
 577         else
 578         {
 579             fepvals->sc_sigma_min = 0;
 580         }
 581     }
 582     if (file_version >= 79)
 583     {
 584         serializer->doBool(&fepvals->bScCoul);
 585     }
 586     else
 587     {
 588         fepvals->bScCoul = TRUE;
 589     }
 590     if (file_version >= 64)
 591     {
 592         serializer->doInt(&fepvals->nstdhdl);
 593     }
 594     else
 595     {
 596         fepvals->nstdhdl = 1;
 597     }
 598
 599     if (file_version >= 73)
 600     {
 601         serializer->doEnumAsInt(&fepvals->separate_dhdl_file);
 602         serializer->doEnumAsInt(&fepvals->dhdl_derivatives);
 603     }
 604     else
 605     {
 606         fepvals->separate_dhdl_file = SeparateDhdlFile::Yes;
 607         fepvals->dhdl_derivatives   = DhDlDerivativeCalculation::Yes;
 608     }
 609     if (file_version >= 71)
 610     {
 611         serializer->doInt(&fepvals->dh_hist_size);
 612         serializer->doDouble(&fepvals->dh_hist_spacing);
 613     }
 614     else
 615     {
 616         fepvals->dh_hist_size    = 0;
 617         fepvals->dh_hist_spacing = 0.1;
 618     }
 619     if (file_version >= 79)
 620     {
 621         serializer->doEnumAsInt(&fepvals->edHdLPrintEnergy);
 622     }
 623     else
 624     {
 625         fepvals->edHdLPrintEnergy = FreeEnergyPrintEnergy::No;
 626     }
 627     if (file_version >= tpxv_SoftcoreGapsys)
 628     {
 629         serializer->doInt(reinterpret_cast<int*>(&fepvals->softcoreFunction));
 630         serializer->doReal(&fepvals->scGapsysScaleLinpointLJ);
 631         serializer->doReal(&fepvals->scGapsysScaleLinpointQ);
 632         serializer->doReal(&fepvals->scGapsysSigmaLJ);
 633     }
 634     else
 635     {
 636         fepvals->softcoreFunction        = SoftcoreType::Beutler;
 637         fepvals->scGapsysScaleLinpointLJ = 0.85;
 638         fepvals->scGapsysScaleLinpointQ  = 0.3;
 639         fepvals->scGapsysSigmaLJ         = 0.3;
 640     }
 641
 642     /* handle lambda_neighbors */
 643     if ((file_version >= 83 && file_version < 90) || file_version >= 92)
 644     {
 645         serializer->doInt(&fepvals->lambda_neighbors);
 646         if ((fepvals->lambda_neighbors >= 0) && (fepvals->init_fep_state >= 0)
 647             && (fepvals->init_lambda < 0))
 648         {
 649             fepvals->lambda_start_n = (fepvals->init_fep_state - fepvals->lambda_neighbors);
 650             fepvals->lambda_stop_n  = (fepvals->init_fep_state + fepvals->lambda_neighbors + 1);
 651             if (fepvals->lambda_start_n < 0)
 652             {
 653                 fepvals->lambda_start_n = 0;
 654             }
 655             if (fepvals->lambda_stop_n >= fepvals->n_lambda)
 656             {
 657                 fepvals->lambda_stop_n = fepvals->n_lambda;
 658             }
 659         }
 660         else
 661         {
 662             fepvals->lambda_start_n = 0;
 663             fepvals->lambda_stop_n  = fepvals->n_lambda;
 664         }
 665     }
 666     else
 667     {
 668         fepvals->lambda_start_n = 0;
 669         fepvals->lambda_stop_n  = fepvals->n_lambda;
 670     }
 671 }
 672
 673 static void do_pull(gmx::ISerializer* serializer, pull_params_t* pull, int file_version, PullingAlgorithm ePullOld)
 674 {
 675     PullGroupGeometry eGeomOld = PullGroupGeometry::Count;
 676     ivec              dimOld;
 677     int               g;
 678
 679     if (file_version >= 95)
 680     {
 681         serializer->doInt(&pull->ngroup);
 682     }
 683     serializer->doInt(&pull->ncoord);
 684     if (file_version < 95)
 685     {
 686         pull->ngroup = pull->ncoord + 1;
 687     }
 688     if (file_version < tpxv_PullCoordTypeGeom)
 689     {
 690         real dum;
 691
 692         serializer->doEnumAsInt(&eGeomOld);
 693         serializer->doIvec(&dimOld);
 694         /* The inner cylinder radius, now removed */
 695         serializer->doReal(&dum);
 696     }
 697     serializer->doReal(&pull->cylinder_r);
 698     serializer->doReal(&pull->constr_tol);
 699     if (file_version >= 95)
 700     {
 701         serializer->doBool(&pull->bPrintCOM);
 702         /* With file_version < 95 this value is set below */
 703     }
 704     if (file_version >= tpxv_ReplacePullPrintCOM12)
 705     {
 706         serializer->doBool(&pull->bPrintRefValue);
 707         serializer->doBool(&pull->bPrintComp);
 708     }
 709     else if (file_version >= tpxv_PullCoordTypeGeom)
 710     {
 711         int idum;
 712         serializer->doInt(&idum); /* used to be bPrintCOM2 */
 713         serializer->doBool(&pull->bPrintRefValue);
 714         serializer->doBool(&pull->bPrintComp);
 715     }
 716     else
 717     {
 718         pull->bPrintRefValue = FALSE;
 719         pull->bPrintComp     = TRUE;
 720     }
 721     serializer->doInt(&pull->nstxout);
 722     serializer->doInt(&pull->nstfout);
 723     if (file_version >= tpxv_PullPrevStepCOMAsReference)
 724     {
 725         serializer->doBool(&pull->bSetPbcRefToPrevStepCOM);
 726     }
 727     else
 728     {
 729         pull->bSetPbcRefToPrevStepCOM = FALSE;
 730     }
 731     pull->group.resize(pull->ngroup);
 732     pull->coord.resize(pull->ncoord);
 733     if (file_version < 95)
 734     {
 735         /* epullgPOS for position pulling, before epullgDIRPBC was removed */
 736         if (eGeomOld == PullGroupGeometry::DirectionPBC)
 737         {
 738             gmx_fatal(FARGS, "pull-geometry=position is no longer supported");
 739         }
 740         if (eGeomOld > PullGroupGeometry::DirectionPBC)
 741         {
 742             switch (eGeomOld)
 743             {
 744                 case (PullGroupGeometry::DirectionRelative):
 745                     eGeomOld = PullGroupGeometry::DirectionPBC;
 746                     break;
 747                 case (PullGroupGeometry::Angle):
 748                     eGeomOld = PullGroupGeometry::DirectionRelative;
 749                     break;
 750                 case (PullGroupGeometry::Dihedral): eGeomOld = PullGroupGeometry::Angle; break;
 751                 case (PullGroupGeometry::AngleAxis): eGeomOld = PullGroupGeometry::Dihedral; break;
 752                 case (PullGroupGeometry::Count): eGeomOld = PullGroupGeometry::AngleAxis; break;
 753                 default: GMX_RELEASE_ASSERT(false, "Unhandled old pull type");
 754             }
 755         }
 756
 757         for (g = 0; g < pull->ngroup; g++)
 758         {
 759             /* We read and ignore a pull coordinate for group 0 */
 760             do_pullgrp_tpx_pre95(serializer, &pull->group[g], &pull->coord[std::max(g - 1, 0)]);
 761             if (g > 0)
 762             {
 763                 pull->coord[g - 1].group[0] = 0;
 764                 pull->coord[g - 1].group[1] = g;
 765             }
 766         }
 767
 768         pull->bPrintCOM = (!pull->group[0].ind.empty());
 769     }
 770     else
 771     {
 772         for (g = 0; g < pull->ngroup; g++)
 773         {
 774             do_pull_group(serializer, &pull->group[g]);
 775         }
 776         for (g = 0; g < pull->ncoord; g++)
 777         {
 778             do_pull_coord(serializer, &pull->coord[g], file_version, ePullOld, eGeomOld, dimOld);
 779             if (serializer->reading())
 780             {
 781                 pull->coord[g].coordIndex = g;
 782             }
 783         }
 784     }
 785     if (file_version >= tpxv_PullAverage)
 786     {
 787         gmx_bool v;
 788
 789         v = pull->bXOutAverage;
 790         serializer->doBool(&v);
 791         pull->bXOutAverage = v;
 792         v                  = pull->bFOutAverage;
 793         serializer->doBool(&v);
 794         pull->bFOutAverage = v;
 795     }
 796 }
 797
 798
 799 static void do_rotgrp(gmx::ISerializer* serializer, t_rotgrp* rotg)
 800 {
 801     serializer->doEnumAsInt(&rotg->eType);
 802     if (serializer->reading())
 803     {
 804         int temp = 0;
 805         serializer->doInt(&temp);
 806         rotg->bMassW = static_cast<bool>(temp);
 807     }
 808     else
 809     {
 810         int temp = static_cast<int>(rotg->bMassW);
 811         serializer->doInt(&temp);
 812     }
 813     serializer->doInt(&rotg->nat);
 814     if (serializer->reading())
 815     {
 816         snew(rotg->ind, rotg->nat);
 817     }
 818     serializer->doIntArray(rotg->ind, rotg->nat);
 819     if (serializer->reading())
 820     {
 821         snew(rotg->x_ref, rotg->nat);
 822     }
 823     serializer->doRvecArray(rotg->x_ref, rotg->nat);
 824     serializer->doRvec(&rotg->inputVec);
 825     serializer->doRvec(&rotg->pivot);
 826     serializer->doReal(&rotg->rate);
 827     serializer->doReal(&rotg->k);
 828     serializer->doReal(&rotg->slab_dist);
 829     serializer->doReal(&rotg->min_gaussian);
 830     serializer->doReal(&rotg->eps);
 831     serializer->doEnumAsInt(&rotg->eFittype);
 832     serializer->doInt(&rotg->PotAngle_nstep);
 833     serializer->doReal(&rotg->PotAngle_step);
 834 }
 835
 836 static void do_rot(gmx::ISerializer* serializer, t_rot* rot)
 837 {
 838     int g;
 839
 840     serializer->doInt(&rot->ngrp);
 841     serializer->doInt(&rot->nstrout);
 842     serializer->doInt(&rot->nstsout);
 843     if (serializer->reading())
 844     {
 845         snew(rot->grp, rot->ngrp);
 846     }
 847     for (g = 0; g < rot->ngrp; g++)
 848     {
 849         do_rotgrp(serializer, &rot->grp[g]);
 850     }
 851 }
 852
 853
 854 static void do_swapgroup(gmx::ISerializer* serializer, t_swapGroup* g)
 855 {
 856
 857     /* Name of the group or molecule */
 858     std::string buf;
 859     if (serializer->reading())
 860     {
 861         serializer->doString(&buf);
 862         g->molname = gmx_strdup(buf.c_str());
 863     }
 864     else
 865     {
 866         buf = g->molname;
 867         serializer->doString(&buf);
 868     }
 869
 870     /* Number of atoms in the group */
 871     serializer->doInt(&g->nat);
 872
 873     /* The group's atom indices */
 874     if (serializer->reading())
 875     {
 876         snew(g->ind, g->nat);
 877     }
 878     serializer->doIntArray(g->ind, g->nat);
 879
 880     /* Requested counts for compartments A and B */
 881     serializer->doIntArray(g->nmolReq.data(), static_cast<int>(Compartment::Count));
 882 }
 883
 884 static void do_swapcoords_tpx(gmx::ISerializer* serializer, t_swapcoords* swap, int file_version)
 885 {
 886     /* Enums for better readability of the code */
 887     enum
 888     {
 889         eCompA = 0,
 890         eCompB
 891     };
 892     enum
 893     {
 894         eChannel0 = 0,
 895         eChannel1
 896     };
 897
 898
 899     if (file_version >= tpxv_CompElPolyatomicIonsAndMultipleIonTypes)
 900     {
 901         /* The total number of swap groups is the sum of the fixed groups
 902          * (split0, split1, solvent), and the user-defined groups (2+ types of ions)
 903          */
 904         serializer->doInt(&swap->ngrp);
 905         if (serializer->reading())
 906         {
 907             snew(swap->grp, swap->ngrp);
 908         }
 909         for (int ig = 0; ig < swap->ngrp; ig++)
 910         {
 911             do_swapgroup(serializer, &swap->grp[ig]);
 912         }
 913         serializer->doBool(&swap->massw_split[eChannel0]);
 914         serializer->doBool(&swap->massw_split[eChannel1]);
 915         serializer->doInt(&swap->nstswap);
 916         serializer->doInt(&swap->nAverage);
 917         serializer->doReal(&swap->threshold);
 918         serializer->doReal(&swap->cyl0r);
 919         serializer->doReal(&swap->cyl0u);
 920         serializer->doReal(&swap->cyl0l);
 921         serializer->doReal(&swap->cyl1r);
 922         serializer->doReal(&swap->cyl1u);
 923         serializer->doReal(&swap->cyl1l);
 924     }
 925     else
 926     {
 927         /*** Support reading older CompEl .tpr files ***/
 928
 929         /* In the original CompEl .tpr files, we always have 5 groups: */
 930         swap->ngrp = 5;
 931         snew(swap->grp, swap->ngrp);
 932
 933         swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].molname = gmx_strdup("split0"); // group 0: split0
 934         swap->grp[static_cast<int>(SwapGroupSplittingType::Split1)].molname = gmx_strdup("split1"); // group 1: split1
 935         swap->grp[static_cast<int>(SwapGroupSplittingType::Solvent)].molname =
 936                 gmx_strdup("solvent");                // group 2: solvent
 937         swap->grp[3].molname = gmx_strdup("anions");  // group 3: anions
 938         swap->grp[4].molname = gmx_strdup("cations"); // group 4: cations
 939
 940         serializer->doInt(&swap->grp[3].nat);
 941         serializer->doInt(&swap->grp[static_cast<int>(SwapGroupSplittingType::Solvent)].nat);
 942         serializer->doInt(&swap->grp[static_cast<int>(SwapGroupSplittingType::Split0)].nat);
 943         serializer->doBool(&swap->massw_split[eChannel0]);
 944         serializer->doInt(&swap->grp[static_cast<int>(SwapGroupSplittingType::Split1)].nat);
 945         serializer->doBool(&swap->massw_split[eChannel1]);
 946         serializer->doInt(&swap->nstswap);
 947         serializer->doInt(&swap->nAverage);
 948         serializer->doReal(&swap->threshold);
 949         serializer->doReal(&swap->cyl0r);
 950         serializer->doReal(&swap->cyl0u);
 951         serializer->doReal(&swap->cyl0l);
 952         serializer->doReal(&swap->cyl1r);
 953         serializer->doReal(&swap->cyl1u);
 954         serializer->doReal(&swap->cyl1l);
 955
 956         // The order[] array keeps compatibility with older .tpr files
 957         // by reading in the groups in the classic order
 958         {
 959             const int order[4] = { 3,
 960                                    static_cast<int>(SwapGroupSplittingType::Solvent),
 961                                    static_cast<int>(SwapGroupSplittingType::Split0),
 962                                    static_cast<int>(SwapGroupSplittingType::Split1) };
 963
 964             for (int ig = 0; ig < 4; ig++)
 965             {
 966                 int g = order[ig];
 967                 snew(swap->grp[g].ind, swap->grp[g].nat);
 968                 serializer->doIntArray(swap->grp[g].ind, swap->grp[g].nat);
 969             }
 970         }
 971
 972         for (int j = eCompA; j <= eCompB; j++)
 973         {
 974             serializer->doInt(&swap->grp[3].nmolReq[j]); // group 3 = anions
 975             serializer->doInt(&swap->grp[4].nmolReq[j]); // group 4 = cations
 976         }
 977     } /* End support reading older CompEl .tpr files */
 978
 979     if (file_version >= tpxv_CompElWithSwapLayerOffset)
 980     {
 981         serializer->doReal(&swap->bulkOffset[eCompA]);
 982         serializer->doReal(&swap->bulkOffset[eCompB]);
 983     }
 984 }
 985
 986 static void do_legacy_efield(gmx::ISerializer* serializer, gmx::KeyValueTreeObjectBuilder* root)
 987 {
 988     const char* const dimName[] = { "x", "y", "z" };
 989
 990     auto appliedForcesObj = root->addObject("applied-forces");
 991     auto efieldObj        = appliedForcesObj.addObject("electric-field");
 992     // The content of the tpr file for this feature has
 993     // been the same since gromacs 4.0 that was used for
 994     // developing.
 995     for (int j = 0; j < DIM; ++j)
 996     {
 997         int n, nt;
 998         serializer->doInt(&n);
 999         serializer->doInt(&nt);
1000         std::vector<real> aa(n + 1), phi(nt + 1), at(nt + 1), phit(nt + 1);
1001         serializer->doRealArray(aa.data(), n);
1002         serializer->doRealArray(phi.data(), n);
1003         serializer->doRealArray(at.data(), nt);
1004         serializer->doRealArray(phit.data(), nt);
1005         if (n > 0)
1006         {
1007             if (n > 1 || nt > 1)
1008             {
1009                 gmx_fatal(FARGS,
1010                           "Can not handle tpr files with more than one electric field term per "
1011                           "direction.");
1012             }
1013             auto dimObj = efieldObj.addObject(dimName[j]);
1014             dimObj.addValue<real>("E0", aa[0]);
1015             dimObj.addValue<real>("omega", at[0]);
1016             dimObj.addValue<real>("t0", phi[0]);
1017             dimObj.addValue<real>("sigma", phit[0]);
1018         }
1019     }
1020 }
1021
1022
1023 static void do_inputrec(gmx::ISerializer* serializer, t_inputrec* ir, int file_version)
1024 {
1025     int      i, j, k, idum = 0;
1026     real     rdum;
1027     gmx_bool bdum = false;
1028
1029     if (file_version != tpx_version)
1030     {
1031         /* Give a warning about features that are not accessible */
1032         fprintf(stderr, "Note: file tpx version %d, software tpx version %d\n", file_version, tpx_version);
1033     }
1034
1035     if (file_version == 0)
1036     {
1037         return;
1038     }
1039
1040     gmx::KeyValueTreeBuilder       paramsBuilder;
1041     gmx::KeyValueTreeObjectBuilder paramsObj = paramsBuilder.rootObject();
1042
1043     /* Basic inputrec stuff */
1044     serializer->doEnumAsInt(&ir->eI);
1045     if (file_version >= 62)
1046     {
1047         serializer->doInt64(&ir->nsteps);
1048     }
1049     else
1050     {
1051         serializer->doInt(&idum);
1052         ir->nsteps = idum;
1053     }
1054
1055     if (file_version >= 62)
1056     {
1057         serializer->doInt64(&ir->init_step);
1058     }
1059     else
1060     {
1061         serializer->doInt(&idum);
1062         ir->init_step = idum;
1063     }
1064
1065     serializer->doInt(&ir->simulation_part);
1066
1067     if (file_version >= tpxv_MTS)
1068     {
1069         serializer->doBool(&ir->useMts);
1070         int numLevels = ir->mtsLevels.size();
1071         if (ir->useMts)
1072         {
1073             serializer->doInt(&numLevels);
1074         }
1075         ir->mtsLevels.resize(numLevels);
1076         for (auto& mtsLevel : ir->mtsLevels)
1077         {
1078             int forceGroups = mtsLevel.forceGroups.to_ulong();
1079             serializer->doInt(&forceGroups);
1080             mtsLevel.forceGroups = std::bitset<static_cast<int>(gmx::MtsForceGroups::Count)>(forceGroups);
1081             serializer->doInt(&mtsLevel.stepFactor);
1082         }
1083     }
1084     else
1085     {
1086         ir->useMts = false;
1087         ir->mtsLevels.clear();
1088     }
1089
1090     if (file_version >= 67)
1091     {
1092         serializer->doInt(&ir->nstcalcenergy);
1093     }
1094     else
1095     {
1096         ir->nstcalcenergy = 1;
1097     }
1098     if (file_version >= 81)
1099     {
1100         serializer->doEnumAsInt(&ir->cutoff_scheme);
1101         if (file_version < 94)
1102         {
1103             // Need to invert the scheme order
1104             switch (ir->cutoff_scheme)
1105             {
1106                 case (CutoffScheme::Group): ir->cutoff_scheme = CutoffScheme::Verlet; break;
1107                 case (CutoffScheme::Verlet): ir->cutoff_scheme = CutoffScheme::Group; break;
1108                 default: GMX_RELEASE_ASSERT(false, "Unhandled cutoff scheme type");
1109             }
1110         }
1111     }
1112     else
1113     {
1114         ir->cutoff_scheme = CutoffScheme::Group;
1115     }
1116     serializer->doInt(&idum); /* used to be ns_type; not used anymore */
1117     serializer->doInt(&ir->nstlist);
1118     serializer->doInt(&idum); /* used to be ndelta; not used anymore */
1119
1120     serializer->doReal(&ir->rtpi);
1121
1122     serializer->doInt(&ir->nstcomm);
1123     serializer->doEnumAsInt(&ir->comm_mode);
1124
1125     /* ignore nstcheckpoint */
1126     if (file_version < tpxv_RemoveObsoleteParameters1)
1127     {
1128         serializer->doInt(&idum);
1129     }
1130
1131     serializer->doInt(&ir->nstcgsteep);
1132
1133     serializer->doInt(&ir->nbfgscorr);
1134
1135     serializer->doInt(&ir->nstlog);
1136     serializer->doInt(&ir->nstxout);
1137     serializer->doInt(&ir->nstvout);
1138     serializer->doInt(&ir->nstfout);
1139     serializer->doInt(&ir->nstenergy);
1140     serializer->doInt(&ir->nstxout_compressed);
1141     if (file_version >= 59)
1142     {
1143         serializer->doDouble(&ir->init_t);
1144         serializer->doDouble(&ir->delta_t);
1145     }
1146     else
1147     {
1148         serializer->doReal(&rdum);
1149         ir->init_t = rdum;
1150         serializer->doReal(&rdum);
1151         ir->delta_t = rdum;
1152     }
1153     serializer->doReal(&ir->x_compression_precision);
1154     if (file_version >= 81)
1155     {
1156         serializer->doReal(&ir->verletbuf_tol);
1157     }
1158     else
1159     {
1160         ir->verletbuf_tol = 0;
1161     }
1162     serializer->doReal(&ir->rlist);
1163     if (file_version >= 67 && file_version < tpxv_RemoveTwinRange)
1164     {
1165         if (serializer->reading())
1166         {
1167             // Reading such a file version could invoke the twin-range
1168             // scheme, about which mdrun should give a fatal error.
1169             real dummy_rlistlong = -1;
1170             serializer->doReal(&dummy_rlistlong);
1171
1172             ir->useTwinRange = (ir->rlist > 0 && (dummy_rlistlong == 0 || dummy_rlistlong > ir->rlist));
1173             // When true, this forces mdrun to issue an error (regardless of
1174             // ir->cutoff_scheme).
1175             //
1176             // Otherwise, grompp used to set rlistlong actively. Users
1177             // were probably also confused and set rlistlong == rlist.
1178             // However, in all remaining cases, it is safe to let
1179             // mdrun proceed normally.
1180         }
1181     }
1182     else
1183     {
1184         // No need to read or write anything
1185         ir->useTwinRange = false;
1186     }
1187     if (file_version >= 82 && file_version != 90)
1188     {
1189         // Multiple time-stepping is no longer enabled, but the old
1190         // support required the twin-range scheme, for which mdrun
1191         // already emits a fatal error.
1192         int dummy_nstcalclr = -1;
1193         serializer->doInt(&dummy_nstcalclr);
1194     }
1195     serializer->doEnumAsInt(&ir->coulombtype);
1196     if (file_version >= 81)
1197     {
1198         serializer->doEnumAsInt(&ir->coulomb_modifier);
1199     }
1200     else
1201     {
1202         ir->coulomb_modifier = (ir->cutoff_scheme == CutoffScheme::Verlet ? InteractionModifiers::PotShift
1203                                                                           : InteractionModifiers::None);
1204     }
1205     serializer->doReal(&ir->rcoulomb_switch);
1206     serializer->doReal(&ir->rcoulomb);
1207     serializer->doEnumAsInt(&ir->vdwtype);
1208     if (file_version >= 81)
1209     {
1210         serializer->doEnumAsInt(&ir->vdw_modifier);
1211     }
1212     else
1213     {
1214         ir->vdw_modifier = (ir->cutoff_scheme == CutoffScheme::Verlet ? InteractionModifiers::PotShift
1215                                                                       : InteractionModifiers::None);
1216     }
1217     serializer->doReal(&ir->rvdw_switch);
1218     serializer->doReal(&ir->rvdw);
1219     serializer->doEnumAsInt(&ir->eDispCorr);
1220     serializer->doReal(&ir->epsilon_r);
1221     serializer->doReal(&ir->epsilon_rf);
1222     serializer->doReal(&ir->tabext);
1223
1224     // This permits reading a .tpr file that used implicit solvent,
1225     // and later permitting mdrun to refuse to run it.
1226     if (serializer->reading())
1227     {
1228         if (file_version < tpxv_RemoveImplicitSolvation)
1229         {
1230             serializer->doInt(&idum);
1231             serializer->doInt(&idum);
1232             serializer->doReal(&rdum);
1233             serializer->doReal(&rdum);
1234             serializer->doInt(&idum);
1235             ir->implicit_solvent = (idum > 0);
1236         }
1237         else
1238         {
1239             ir->implicit_solvent = false;
1240         }
1241         if (file_version < tpxv_RemoveImplicitSolvation)
1242         {
1243             serializer->doReal(&rdum);
1244             serializer->doReal(&rdum);
1245             serializer->doReal(&rdum);
1246             serializer->doReal(&rdum);
1247             if (file_version >= 60)
1248             {
1249                 serializer->doReal(&rdum);
1250                 serializer->doInt(&idum);
1251             }
1252             serializer->doReal(&rdum);
1253         }
1254     }
1255
1256     if (file_version >= 81)
1257     {
1258         serializer->doReal(&ir->fourier_spacing);
1259     }
1260     else
1261     {
1262         ir->fourier_spacing = 0.0;
1263     }
1264     serializer->doInt(&ir->nkx);
1265     serializer->doInt(&ir->nky);
1266     serializer->doInt(&ir->nkz);
1267     serializer->doInt(&ir->pme_order);
1268     serializer->doReal(&ir->ewald_rtol);
1269
1270     if (file_version >= 93)
1271     {
1272         serializer->doReal(&ir->ewald_rtol_lj);
1273     }
1274     else
1275     {
1276         ir->ewald_rtol_lj = ir->ewald_rtol;
1277     }
1278     serializer->doEnumAsInt(&ir->ewald_geometry);
1279     serializer->doReal(&ir->epsilon_surface);
1280
1281     /* ignore bOptFFT */
1282     if (file_version < tpxv_RemoveObsoleteParameters1)
1283     {
1284         serializer->doBool(&bdum);
1285     }
1286
1287     if (file_version >= 93)
1288     {
1289         serializer->doEnumAsInt(&ir->ljpme_combination_rule);
1290     }
1291     serializer->doBool(&ir->bContinuation);
1292     serializer->doEnumAsInt(&ir->etc);
1293     /* before version 18, ir->etc was a gmx_bool (ir->btc),
1294      * but the values 0 and 1 still mean no and
1295      * berendsen temperature coupling, respectively.
1296      */
1297     if (file_version >= 79)
1298     {
1299         serializer->doBool(&ir->bPrintNHChains);
1300     }
1301     if (file_version >= 71)
1302     {
1303         serializer->doInt(&ir->nsttcouple);
1304     }
1305     else
1306     {
1307         ir->nsttcouple = ir->nstcalcenergy;
1308     }
1309     serializer->doEnumAsInt(&ir->epc);
1310     serializer->doEnumAsInt(&ir->epct);
1311     if (file_version >= 71)
1312     {
1313         serializer->doInt(&ir->nstpcouple);
1314     }
1315     else
1316     {
1317         ir->nstpcouple = ir->nstcalcenergy;
1318     }
1319     serializer->doReal(&ir->tau_p);
1320     serializer->doRvec(&ir->ref_p[XX]);
1321     serializer->doRvec(&ir->ref_p[YY]);
1322     serializer->doRvec(&ir->ref_p[ZZ]);
1323     serializer->doRvec(&ir->compress[XX]);
1324     serializer->doRvec(&ir->compress[YY]);
1325     serializer->doRvec(&ir->compress[ZZ]);
1326     serializer->doEnumAsInt(&ir->refcoord_scaling);
1327     serializer->doRvec(&ir->posres_com);
1328     serializer->doRvec(&ir->posres_comB);
1329
1330     if (file_version < 79)
1331     {
1332         serializer->doInt(&ir->andersen_seed);
1333     }
1334     else
1335     {
1336         ir->andersen_seed = 0;
1337     }
1338
1339     serializer->doReal(&ir->shake_tol);
1340
1341     serializer->doEnumAsInt(&ir->efep);
1342     if (serializer->reading())
1343     {
1344         ir->fepvals = std::make_unique<t_lambda>();
1345     }
1346     do_fepvals(serializer, ir->fepvals.get(), file_version);
1347
1348     if (file_version >= 79)
1349     {
1350         serializer->doBool(&ir->bSimTemp);
1351         if (ir->bSimTemp)
1352         {
1353             ir->bSimTemp = TRUE;
1354         }
1355     }
1356     else
1357     {
1358         ir->bSimTemp = FALSE;
1359     }
1360     if (ir->bSimTemp)
1361     {
1362         if (serializer->reading())
1363         {
1364             ir->simtempvals = std::make_unique<t_simtemp>();
1365         }
1366         do_simtempvals(serializer, ir->simtempvals.get(), ir->fepvals->n_lambda, file_version);
1367     }
1368
1369     if (file_version >= 79)
1370     {
1371         serializer->doBool(&ir->bExpanded);
1372     }
1373     if (ir->bExpanded)
1374     {
1375         if (serializer->reading())
1376         {
1377             ir->expandedvals = std::make_unique<t_expanded>();
1378         }
1379         do_expandedvals(serializer, ir->expandedvals.get(), ir->fepvals.get(), file_version);
1380     }
1381
1382     serializer->doEnumAsInt(&ir->eDisre);
1383     serializer->doEnumAsInt(&ir->eDisreWeighting);
1384     serializer->doBool(&ir->bDisreMixed);
1385     serializer->doReal(&ir->dr_fc);
1386     serializer->doReal(&ir->dr_tau);
1387     serializer->doInt(&ir->nstdisreout);
1388     serializer->doReal(&ir->orires_fc);
1389     serializer->doReal(&ir->orires_tau);
1390     serializer->doInt(&ir->nstorireout);
1391
1392     /* ignore dihre_fc */
1393     if (file_version < 79)
1394     {
1395         serializer->doReal(&rdum);
1396     }
1397
1398     serializer->doReal(&ir->em_stepsize);
1399     serializer->doReal(&ir->em_tol);
1400     serializer->doBool(&ir->bShakeSOR);
1401     serializer->doInt(&ir->niter);
1402     serializer->doReal(&ir->fc_stepsize);
1403     serializer->doEnumAsInt(&ir->eConstrAlg);
1404     serializer->doInt(&ir->nProjOrder);
1405     serializer->doReal(&ir->LincsWarnAngle);
1406     serializer->doInt(&ir->nLincsIter);
1407     serializer->doReal(&ir->bd_fric);
1408     if (file_version >= tpxv_Use64BitRandomSeed)
1409     {
1410         serializer->doInt64(&ir->ld_seed);
1411     }
1412     else
1413     {
1414         serializer->doInt(&idum);
1415         ir->ld_seed = idum;
1416     }
1417
1418     for (i = 0; i < DIM; i++)
1419     {
1420         serializer->doRvec(&ir->deform[i]);
1421     }
1422     serializer->doReal(&ir->cos_accel);
1423
1424     serializer->doInt(&ir->userint1);
1425     serializer->doInt(&ir->userint2);
1426     serializer->doInt(&ir->userint3);
1427     serializer->doInt(&ir->userint4);
1428     serializer->doReal(&ir->userreal1);
1429     serializer->doReal(&ir->userreal2);
1430     serializer->doReal(&ir->userreal3);
1431     serializer->doReal(&ir->userreal4);
1432
1433     /* AdResS is removed, but we need to be able to read old files,
1434        and let mdrun refuse to run them */
1435     if (file_version >= 77 && file_version < tpxv_RemoveAdress)
1436     {
1437         serializer->doBool(&ir->bAdress);
1438         if (ir->bAdress)
1439         {
1440             int  idum, numThermoForceGroups, numEnergyGroups;
1441             real rdum;
1442             rvec rvecdum;
1443             serializer->doInt(&idum);
1444             serializer->doReal(&rdum);
1445             serializer->doReal(&rdum);
1446             serializer->doReal(&rdum);
1447             serializer->doInt(&idum);
1448             serializer->doInt(&idum);
1449             serializer->doRvec(&rvecdum);
1450             serializer->doInt(&numThermoForceGroups);
1451             serializer->doReal(&rdum);
1452             serializer->doInt(&numEnergyGroups);
1453             serializer->doInt(&idum);
1454
1455             if (numThermoForceGroups > 0)
1456             {
1457                 std::vector<int> idumn(numThermoForceGroups);
1458                 serializer->doIntArray(idumn.data(), idumn.size());
1459             }
1460             if (numEnergyGroups > 0)
1461             {
1462                 std::vector<int> idumn(numEnergyGroups);
1463                 serializer->doIntArray(idumn.data(), idumn.size());
1464             }
1465         }
1466     }
1467     else
1468     {
1469         ir->bAdress = FALSE;
1470     }
1471
1472     /* pull stuff */
1473     {
1474         PullingAlgorithm ePullOld = PullingAlgorithm::Umbrella;
1475
1476         if (file_version >= tpxv_PullCoordTypeGeom)
1477         {
1478             serializer->doBool(&ir->bPull);
1479         }
1480         else
1481         {
1482             serializer->doEnumAsInt(&ePullOld);
1483             ir->bPull = (ePullOld != PullingAlgorithm::Umbrella);
1484             /* We removed the first ePull=ePullNo for the enum */
1485             switch (ePullOld)
1486             {
1487                 case (PullingAlgorithm::Umbrella): break; // this is equal to not using pulling
1488                 case (PullingAlgorithm::Constraint): ePullOld = PullingAlgorithm::Umbrella; break;
1489                 case (PullingAlgorithm::ConstantForce):
1490                     ePullOld = PullingAlgorithm::Constraint;
1491                     break;
1492                 case (PullingAlgorithm::FlatBottom):
1493                     ePullOld = PullingAlgorithm::ConstantForce;
1494                     break;
1495                 case (PullingAlgorithm::FlatBottomHigh):
1496                     ePullOld = PullingAlgorithm::FlatBottom;
1497                     break;
1498                 case (PullingAlgorithm::External):
1499                     ePullOld = PullingAlgorithm::FlatBottomHigh;
1500                     break;
1501                 case (PullingAlgorithm::Count): ePullOld = PullingAlgorithm::External; break;
1502                 default: GMX_RELEASE_ASSERT(false, "Unhandled old pull algorithm");
1503             }
1504         }
1505         if (ir->bPull)
1506         {
1507             if (serializer->reading())
1508             {
1509                 ir->pull = std::make_unique<pull_params_t>();
1510             }
1511             do_pull(serializer, ir->pull.get(), file_version, ePullOld);
1512         }
1513     }
1514
1515     if (file_version >= tpxv_AcceleratedWeightHistogram)
1516     {
1517         serializer->doBool(&ir->bDoAwh);
1518
1519         if (ir->bDoAwh)
1520         {
1521             if (serializer->reading())
1522             {
1523                 ir->awhParams = std::make_unique<gmx::AwhParams>(serializer);
1524             }
1525             else
1526             {
1527                 ir->awhParams->serialize(serializer);
1528             }
1529         }
1530     }
1531     else
1532     {
1533         ir->bDoAwh = FALSE;
1534     }
1535
1536     /* Enforced rotation */
1537     if (file_version >= 74)
1538     {
1539         serializer->doBool(&ir->bRot);
1540         if (ir->bRot)
1541         {
1542             if (serializer->reading())
1543             {
1544                 snew(ir->rot, 1);
1545             }
1546             do_rot(serializer, ir->rot);
1547         }
1548     }
1549     else
1550     {
1551         ir->bRot = FALSE;
1552     }
1553
1554     /* Interactive molecular dynamics */
1555     if (file_version >= tpxv_InteractiveMolecularDynamics)
1556     {
1557         serializer->doBool(&ir->bIMD);
1558         if (ir->bIMD)
1559         {
1560             if (serializer->reading())
1561             {
1562                 snew(ir->imd, 1);
1563             }
1564             do_imd(serializer, ir->imd);
1565         }
1566     }
1567     else
1568     {
1569         /* We don't support IMD sessions for old .tpr files */
1570         ir->bIMD = FALSE;
1571     }
1572
1573     /* grpopts stuff */
1574     serializer->doInt(&ir->opts.ngtc);
1575     if (file_version >= 69)
1576     {
1577         serializer->doInt(&ir->opts.nhchainlength);
1578     }
1579     else
1580     {
1581         ir->opts.nhchainlength = 1;
1582     }
1583     if (serializer->reading() && file_version >= tpxv_RemovedConstantAcceleration
1584         && file_version < tpxv_ReaddedConstantAcceleration)
1585     {
1586         ir->opts.ngacc = 0;
1587     }
1588     else
1589     {
1590         serializer->doInt(&ir->opts.ngacc);
1591     }
1592     serializer->doInt(&ir->opts.ngfrz);
1593     serializer->doInt(&ir->opts.ngener);
1594
1595     if (serializer->reading())
1596     {
1597         snew(ir->opts.nrdf, ir->opts.ngtc);
1598         snew(ir->opts.ref_t, ir->opts.ngtc);
1599         snew(ir->opts.annealing, ir->opts.ngtc);
1600         snew(ir->opts.anneal_npoints, ir->opts.ngtc);
1601         snew(ir->opts.anneal_time, ir->opts.ngtc);
1602         snew(ir->opts.anneal_temp, ir->opts.ngtc);
1603         snew(ir->opts.tau_t, ir->opts.ngtc);
1604         snew(ir->opts.nFreeze, ir->opts.ngfrz);
1605         snew(ir->opts.acceleration, ir->opts.ngacc);
1606         snew(ir->opts.egp_flags, ir->opts.ngener * ir->opts.ngener);
1607     }
1608     if (ir->opts.ngtc > 0)
1609     {
1610         serializer->doRealArray(ir->opts.nrdf, ir->opts.ngtc);
1611         serializer->doRealArray(ir->opts.ref_t, ir->opts.ngtc);
1612         serializer->doRealArray(ir->opts.tau_t, ir->opts.ngtc);
1613     }
1614     if (ir->opts.ngfrz > 0)
1615     {
1616         serializer->doIvecArray(ir->opts.nFreeze, ir->opts.ngfrz);
1617     }
1618     if (ir->opts.ngacc > 0)
1619     {
1620         serializer->doRvecArray(ir->opts.acceleration, ir->opts.ngacc);
1621     }
1622     if (serializer->reading())
1623     {
1624         ir->useConstantAcceleration = false;
1625         for (int g = 0; g < ir->opts.ngacc; g++)
1626         {
1627             if (norm2(ir->opts.acceleration[g]) != 0)
1628             {
1629                 ir->useConstantAcceleration = true;
1630             }
1631         }
1632     }
1633     serializer->doIntArray(ir->opts.egp_flags, ir->opts.ngener * ir->opts.ngener);
1634
1635     /* First read the lists with annealing and npoints for each group */
1636     serializer->doEnumArrayAsInt(ir->opts.annealing, ir->opts.ngtc);
1637     serializer->doIntArray(ir->opts.anneal_npoints, ir->opts.ngtc);
1638     for (j = 0; j < (ir->opts.ngtc); j++)
1639     {
1640         k = ir->opts.anneal_npoints[j];
1641         if (serializer->reading())
1642         {
1643             snew(ir->opts.anneal_time[j], k);
1644             snew(ir->opts.anneal_temp[j], k);
1645         }
1646         serializer->doRealArray(ir->opts.anneal_time[j], k);
1647         serializer->doRealArray(ir->opts.anneal_temp[j], k);
1648     }
1649     /* Walls */
1650     {
1651         serializer->doInt(&ir->nwall);
1652         serializer->doEnumAsInt(&ir->wall_type);
1653         serializer->doReal(&ir->wall_r_linpot);
1654         serializer->doInt(&ir->wall_atomtype[0]);
1655         serializer->doInt(&ir->wall_atomtype[1]);
1656         serializer->doReal(&ir->wall_density[0]);
1657         serializer->doReal(&ir->wall_density[1]);
1658         serializer->doReal(&ir->wall_ewald_zfac);
1659     }
1660
1661     /* Cosine stuff for electric fields */
1662     if (file_version < tpxv_GenericParamsForElectricField)
1663     {
1664         do_legacy_efield(serializer, &paramsObj);
1665     }
1666
1667     /* Swap ions */
1668     if (file_version >= tpxv_ComputationalElectrophysiology)
1669     {
1670         serializer->doEnumAsInt(&ir->eSwapCoords);
1671         if (ir->eSwapCoords != SwapType::No)
1672         {
1673             if (serializer->reading())
1674             {
1675                 snew(ir->swap, 1);
1676             }
1677             do_swapcoords_tpx(serializer, ir->swap, file_version);
1678         }
1679     }
1680
1681     /* QMMM reading - despite defunct we require reading for backwards
1682      * compability and correct serialization
1683      */
1684     {
1685         serializer->doBool(&ir->bQMMM);
1686         int qmmmScheme;
1687         serializer->doInt(&qmmmScheme);
1688         real unusedScalefactor;
1689         serializer->doReal(&unusedScalefactor);
1690
1691         // this is still used in Mimic
1692         serializer->doInt(&ir->opts.ngQM);
1693         if (ir->opts.ngQM > 0 && ir->bQMMM)
1694         {
1695             /* We leave in dummy i/o for removed parameters to avoid
1696              * changing the tpr format.
1697              */
1698             std::vector<int> dummyIntVec(4 * ir->opts.ngQM, 0);
1699             serializer->doIntArray(dummyIntVec.data(), dummyIntVec.size());
1700             dummyIntVec.clear();
1701
1702             // std::vector<bool> has no data()
1703             std::vector<char> dummyBoolVec(ir->opts.ngQM * sizeof(bool) / sizeof(char));
1704             serializer->doBoolArray(reinterpret_cast<bool*>(dummyBoolVec.data()), dummyBoolVec.size());
1705             dummyBoolVec.clear();
1706
1707             dummyIntVec.resize(2 * ir->opts.ngQM, 0);
1708             serializer->doIntArray(dummyIntVec.data(), dummyIntVec.size());
1709             dummyIntVec.clear();
1710
1711             std::vector<real> dummyRealVec(2 * ir->opts.ngQM, 0);
1712             serializer->doRealArray(dummyRealVec.data(), dummyRealVec.size());
1713             dummyRealVec.clear();
1714
1715             dummyIntVec.resize(3 * ir->opts.ngQM, 0);
1716             serializer->doIntArray(dummyIntVec.data(), dummyIntVec.size());
1717             dummyIntVec.clear();
1718         }
1719         /* end of QMMM stuff */
1720     }
1721
1722     if (file_version >= tpxv_GenericParamsForElectricField)
1723     {
1724         if (serializer->reading())
1725         {
1726             paramsObj.mergeObject(gmx::deserializeKeyValueTree(serializer));
1727         }
1728         else
1729         {
1730             GMX_RELEASE_ASSERT(ir->params != nullptr,
1731                                "Parameters should be present when writing inputrec");
1732             gmx::serializeKeyValueTree(*ir->params, serializer);
1733         }
1734     }
1735     if (serializer->reading())
1736     {
1737         ir->params = new gmx::KeyValueTreeObject(paramsBuilder.build());
1738         // Initialize internal parameters to an empty kvt for all tpr versions
1739         ir->internalParameters = std::make_unique<gmx::KeyValueTreeObject>();
1740     }
1741
1742     if (file_version >= tpxv_GenericInternalParameters)
1743     {
1744         if (serializer->reading())
1745         {
1746             ir->internalParameters =
1747                     std::make_unique<gmx::KeyValueTreeObject>(gmx::deserializeKeyValueTree(serializer));
1748         }
1749         else
1750         {
1751             GMX_RELEASE_ASSERT(ir->internalParameters != nullptr,
1752                                "Parameters should be present when writing inputrec");
1753             gmx::serializeKeyValueTree(*ir->internalParameters, serializer);
1754         }
1755     }
1756 }
1757
1758
1759 static void do_harm(gmx::ISerializer* serializer, t_iparams* iparams)
1760 {
1761     serializer->doReal(&iparams->harmonic.rA);
1762     serializer->doReal(&iparams->harmonic.krA);
1763     serializer->doReal(&iparams->harmonic.rB);
1764     serializer->doReal(&iparams->harmonic.krB);
1765 }
1766
1767 static void do_iparams(gmx::ISerializer* serializer, t_functype ftype, t_iparams* iparams, int file_version)
1768 {
1769     int  idum;
1770     real rdum;
1771
1772     switch (ftype)
1773     {
1774         case F_ANGLES:
1775         case F_G96ANGLES:
1776         case F_BONDS:
1777         case F_G96BONDS:
1778         case F_HARMONIC:
1779         case F_IDIHS:
1780             do_harm(serializer, iparams);
1781             if ((ftype == F_ANGRES || ftype == F_ANGRESZ) && serializer->reading())
1782             {
1783                 /* Correct incorrect storage of parameters */
1784                 iparams->pdihs.phiB = iparams->pdihs.phiA;
1785                 iparams->pdihs.cpB  = iparams->pdihs.cpA;
1786             }
1787             break;
1788         case F_RESTRANGLES:
1789             serializer->doReal(&iparams->harmonic.rA);
1790             serializer->doReal(&iparams->harmonic.krA);
1791             break;
1792         case F_LINEAR_ANGLES:
1793             serializer->doReal(&iparams->linangle.klinA);
1794             serializer->doReal(&iparams->linangle.aA);
1795             serializer->doReal(&iparams->linangle.klinB);
1796             serializer->doReal(&iparams->linangle.aB);
1797             break;
1798         case F_FENEBONDS:
1799             serializer->doReal(&iparams->fene.bm);
1800             serializer->doReal(&iparams->fene.kb);
1801             break;
1802
1803         case F_RESTRBONDS:
1804             serializer->doReal(&iparams->restraint.lowA);
1805             serializer->doReal(&iparams->restraint.up1A);
1806             serializer->doReal(&iparams->restraint.up2A);
1807             serializer->doReal(&iparams->restraint.kA);
1808             serializer->doReal(&iparams->restraint.lowB);
1809             serializer->doReal(&iparams->restraint.up1B);
1810             serializer->doReal(&iparams->restraint.up2B);
1811             serializer->doReal(&iparams->restraint.kB);
1812             break;
1813         case F_TABBONDS:
1814         case F_TABBONDSNC:
1815         case F_TABANGLES:
1816         case F_TABDIHS:
1817             serializer->doReal(&iparams->tab.kA);
1818             serializer->doInt(&iparams->tab.table);
1819             serializer->doReal(&iparams->tab.kB);
1820             break;
1821         case F_CROSS_BOND_BONDS:
1822             serializer->doReal(&iparams->cross_bb.r1e);
1823             serializer->doReal(&iparams->cross_bb.r2e);
1824             serializer->doReal(&iparams->cross_bb.krr);
1825             break;
1826         case F_CROSS_BOND_ANGLES:
1827             serializer->doReal(&iparams->cross_ba.r1e);
1828             serializer->doReal(&iparams->cross_ba.r2e);
1829             serializer->doReal(&iparams->cross_ba.r3e);
1830             serializer->doReal(&iparams->cross_ba.krt);
1831             break;
1832         case F_UREY_BRADLEY:
1833             serializer->doReal(&iparams->u_b.thetaA);
1834             serializer->doReal(&iparams->u_b.kthetaA);
1835             serializer->doReal(&iparams->u_b.r13A);
1836             serializer->doReal(&iparams->u_b.kUBA);
1837             if (file_version >= 79)
1838             {
1839                 serializer->doReal(&iparams->u_b.thetaB);
1840                 serializer->doReal(&iparams->u_b.kthetaB);
1841                 serializer->doReal(&iparams->u_b.r13B);
1842                 serializer->doReal(&iparams->u_b.kUBB);
1843             }
1844             else
1845             {
1846                 iparams->u_b.thetaB  = iparams->u_b.thetaA;
1847                 iparams->u_b.kthetaB = iparams->u_b.kthetaA;
1848                 iparams->u_b.r13B    = iparams->u_b.r13A;
1849                 iparams->u_b.kUBB    = iparams->u_b.kUBA;
1850             }
1851             break;
1852         case F_QUARTIC_ANGLES:
1853             serializer->doReal(&iparams->qangle.theta);
1854             serializer->doRealArray(iparams->qangle.c, 5);
1855             break;
1856         case F_BHAM:
1857             serializer->doReal(&iparams->bham.a);
1858             serializer->doReal(&iparams->bham.b);
1859             serializer->doReal(&iparams->bham.c);
1860             break;
1861         case F_MORSE:
1862             serializer->doReal(&iparams->morse.b0A);
1863             serializer->doReal(&iparams->morse.cbA);
1864             serializer->doReal(&iparams->morse.betaA);
1865             if (file_version >= 79)
1866             {
1867                 serializer->doReal(&iparams->morse.b0B);
1868                 serializer->doReal(&iparams->morse.cbB);
1869                 serializer->doReal(&iparams->morse.betaB);
1870             }
1871             else
1872             {
1873                 iparams->morse.b0B   = iparams->morse.b0A;
1874                 iparams->morse.cbB   = iparams->morse.cbA;
1875                 iparams->morse.betaB = iparams->morse.betaA;
1876             }
1877             break;
1878         case F_CUBICBONDS:
1879             serializer->doReal(&iparams->cubic.b0);
1880             serializer->doReal(&iparams->cubic.kb);
1881             serializer->doReal(&iparams->cubic.kcub);
1882             break;
1883         case F_CONNBONDS: break;
1884         case F_POLARIZATION: serializer->doReal(&iparams->polarize.alpha); break;
1885         case F_ANHARM_POL:
1886             serializer->doReal(&iparams->anharm_polarize.alpha);
1887             serializer->doReal(&iparams->anharm_polarize.drcut);
1888             serializer->doReal(&iparams->anharm_polarize.khyp);
1889             break;
1890         case F_WATER_POL:
1891             serializer->doReal(&iparams->wpol.al_x);
1892             serializer->doReal(&iparams->wpol.al_y);
1893             serializer->doReal(&iparams->wpol.al_z);
1894             serializer->doReal(&iparams->wpol.rOH);
1895             serializer->doReal(&iparams->wpol.rHH);
1896             serializer->doReal(&iparams->wpol.rOD);
1897             break;
1898         case F_THOLE_POL:
1899             serializer->doReal(&iparams->thole.a);
1900             serializer->doReal(&iparams->thole.alpha1);
1901             serializer->doReal(&iparams->thole.alpha2);
1902             serializer->doReal(&iparams->thole.rfac);
1903             break;
1904         case F_LJ:
1905             serializer->doReal(&iparams->lj.c6);
1906             serializer->doReal(&iparams->lj.c12);
1907             break;
1908         case F_LJ14:
1909             serializer->doReal(&iparams->lj14.c6A);
1910             serializer->doReal(&iparams->lj14.c12A);
1911             serializer->doReal(&iparams->lj14.c6B);
1912             serializer->doReal(&iparams->lj14.c12B);
1913             break;
1914         case F_LJC14_Q:
1915             serializer->doReal(&iparams->ljc14.fqq);
1916             serializer->doReal(&iparams->ljc14.qi);
1917             serializer->doReal(&iparams->ljc14.qj);
1918             serializer->doReal(&iparams->ljc14.c6);
1919             serializer->doReal(&iparams->ljc14.c12);
1920             break;
1921         case F_LJC_PAIRS_NB:
1922             serializer->doReal(&iparams->ljcnb.qi);
1923             serializer->doReal(&iparams->ljcnb.qj);
1924             serializer->doReal(&iparams->ljcnb.c6);
1925             serializer->doReal(&iparams->ljcnb.c12);
1926             break;
1927         case F_PDIHS:
1928         case F_PIDIHS:
1929         case F_ANGRES:
1930         case F_ANGRESZ:
1931             serializer->doReal(&iparams->pdihs.phiA);
1932             serializer->doReal(&iparams->pdihs.cpA);
1933             serializer->doReal(&iparams->pdihs.phiB);
1934             serializer->doReal(&iparams->pdihs.cpB);
1935             serializer->doInt(&iparams->pdihs.mult);
1936             break;
1937         case F_RESTRDIHS:
1938             serializer->doReal(&iparams->pdihs.phiA);
1939             serializer->doReal(&iparams->pdihs.cpA);
1940             break;
1941         case F_DISRES:
1942             serializer->doInt(&iparams->disres.label);
1943             serializer->doInt(&iparams->disres.type);
1944             serializer->doReal(&iparams->disres.low);
1945             serializer->doReal(&iparams->disres.up1);
1946             serializer->doReal(&iparams->disres.up2);
1947             serializer->doReal(&iparams->disres.kfac);
1948             break;
1949         case F_ORIRES:
1950             serializer->doInt(&iparams->orires.ex);
1951             serializer->doInt(&iparams->orires.label);
1952             serializer->doInt(&iparams->orires.power);
1953             serializer->doReal(&iparams->orires.c);
1954             serializer->doReal(&iparams->orires.obs);
1955             serializer->doReal(&iparams->orires.kfac);
1956             break;
1957         case F_DIHRES:
1958             if (file_version < 82)
1959             {
1960                 serializer->doInt(&idum);
1961                 serializer->doInt(&idum);
1962             }
1963             serializer->doReal(&iparams->dihres.phiA);
1964             serializer->doReal(&iparams->dihres.dphiA);
1965             serializer->doReal(&iparams->dihres.kfacA);
1966             if (file_version >= 82)
1967             {
1968                 serializer->doReal(&iparams->dihres.phiB);
1969                 serializer->doReal(&iparams->dihres.dphiB);
1970                 serializer->doReal(&iparams->dihres.kfacB);
1971             }
1972             else
1973             {
1974                 iparams->dihres.phiB  = iparams->dihres.phiA;
1975                 iparams->dihres.dphiB = iparams->dihres.dphiA;
1976                 iparams->dihres.kfacB = iparams->dihres.kfacA;
1977             }
1978             break;
1979         case F_POSRES:
1980             serializer->doRvec(&iparams->posres.pos0A);
1981             serializer->doRvec(&iparams->posres.fcA);
1982             serializer->doRvec(&iparams->posres.pos0B);
1983             serializer->doRvec(&iparams->posres.fcB);
1984             break;
1985         case F_FBPOSRES:
1986             serializer->doInt(&iparams->fbposres.geom);
1987             serializer->doRvec(&iparams->fbposres.pos0);
1988             serializer->doReal(&iparams->fbposres.r);
1989             serializer->doReal(&iparams->fbposres.k);
1990             break;
1991         case F_CBTDIHS: serializer->doRealArray(iparams->cbtdihs.cbtcA, NR_CBTDIHS); break;
1992         case F_RBDIHS:
1993             // Fall-through intended
1994         case F_FOURDIHS:
1995             /* Fourier dihedrals are internally represented
1996              * as Ryckaert-Bellemans since those are faster to compute.
1997              */
1998             serializer->doRealArray(iparams->rbdihs.rbcA, NR_RBDIHS);
1999             serializer->doRealArray(iparams->rbdihs.rbcB, NR_RBDIHS);
2000             break;
2001         case F_CONSTR:
2002         case F_CONSTRNC:
2003             serializer->doReal(&iparams->constr.dA);
2004             serializer->doReal(&iparams->constr.dB);
2005             break;
2006         case F_SETTLE:
2007             serializer->doReal(&iparams->settle.doh);
2008             serializer->doReal(&iparams->settle.dhh);
2009             break;
2010         case F_VSITE1: break; // VSite1 has 0 parameters
2011         case F_VSITE2:
2012         case F_VSITE2FD: serializer->doReal(&iparams->vsite.a); break;
2013         case F_VSITE3:
2014         case F_VSITE3FD:
2015         case F_VSITE3FAD:
2016             serializer->doReal(&iparams->vsite.a);
2017             serializer->doReal(&iparams->vsite.b);
2018             break;
2019         case F_VSITE3OUT:
2020         case F_VSITE4FD:
2021         case F_VSITE4FDN:
2022             serializer->doReal(&iparams->vsite.a);
2023             serializer->doReal(&iparams->vsite.b);
2024             serializer->doReal(&iparams->vsite.c);
2025             break;
2026         case F_VSITEN:
2027             serializer->doInt(&iparams->vsiten.n);
2028             serializer->doReal(&iparams->vsiten.a);
2029             break;
2030         case F_GB12_NOLONGERUSED:
2031         case F_GB13_NOLONGERUSED:
2032         case F_GB14_NOLONGERUSED:
2033             // Implicit solvent parameters can still be read, but never used
2034             if (serializer->reading())
2035             {
2036                 if (file_version < 68)
2037                 {
2038                     serializer->doReal(&rdum);
2039                     serializer->doReal(&rdum);
2040                     serializer->doReal(&rdum);
2041                     serializer->doReal(&rdum);
2042                 }
2043                 if (file_version < tpxv_RemoveImplicitSolvation)
2044                 {
2045                     serializer->doReal(&rdum);
2046                     serializer->doReal(&rdum);
2047                     serializer->doReal(&rdum);
2048                     serializer->doReal(&rdum);
2049                     serializer->doReal(&rdum);
2050                 }
2051             }
2052             break;
2053         case F_CMAP:
2054             serializer->doInt(&iparams->cmap.cmapA);
2055             serializer->doInt(&iparams->cmap.cmapB);
2056             break;
2057         default:
2058             gmx_fatal(FARGS,
2059                       "unknown function type %d (%s) in %s line %d",
2060                       ftype,
2061                       interaction_function[ftype].name,
2062                       __FILE__,
2063                       __LINE__);
2064     }
2065 }
2066
2067 static void do_ilist(gmx::ISerializer* serializer, InteractionList* ilist)
2068 {
2069     int nr = ilist->size();
2070     serializer->doInt(&nr);
2071     if (serializer->reading())
2072     {
2073         ilist->iatoms.resize(nr);
2074     }
2075     serializer->doIntArray(ilist->iatoms.data(), ilist->size());
2076 }
2077
2078 static void do_ffparams(gmx::ISerializer* serializer, gmx_ffparams_t* ffparams, int file_version)
2079 {
2080     serializer->doInt(&ffparams->atnr);
2081     int numTypes = ffparams->numTypes();
2082     serializer->doInt(&numTypes);
2083     if (serializer->reading())
2084     {
2085         ffparams->functype.resize(numTypes);
2086         ffparams->iparams.resize(numTypes);
2087     }
2088     /* Read/write all the function types */
2089     serializer->doIntArray(ffparams->functype.data(), ffparams->functype.size());
2090
2091     if (file_version >= 66)
2092     {
2093         serializer->doDouble(&ffparams->reppow);
2094     }
2095     else
2096     {
2097         ffparams->reppow = 12.0;
2098     }
2099
2100     serializer->doReal(&ffparams->fudgeQQ);
2101
2102     /* Check whether all these function types are supported by the code.
2103      * In practice the code is backwards compatible, which means that the
2104      * numbering may have to be altered from old numbering to new numbering
2105      */
2106     for (int i = 0; i < ffparams->numTypes(); i++)
2107     {
2108         if (serializer->reading())
2109         {
2110             /* Loop over file versions */
2111             for (int k = 0; k < NFTUPD; k++)
2112             {
2113                 /* Compare the read file_version to the update table */
2114                 if ((file_version < ftupd[k].fvnr) && (ffparams->functype[i] >= ftupd[k].ftype))
2115                 {
2116                     ffparams->functype[i] += 1;
2117                 }
2118             }
2119         }
2120
2121         do_iparams(serializer, ffparams->functype[i], &ffparams->iparams[i], file_version);
2122     }
2123 }
2124
2125 static void add_settle_atoms(InteractionList* ilist)
2126 {
2127     int i;
2128
2129     /* Settle used to only store the first atom: add the other two */
2130     ilist->iatoms.resize(2 * ilist->size());
2131     for (i = ilist->size() / 4 - 1; i >= 0; i--)
2132     {
2133         ilist->iatoms[4 * i + 0] = ilist->iatoms[2 * i + 0];
2134         ilist->iatoms[4 * i + 1] = ilist->iatoms[2 * i + 1];
2135         ilist->iatoms[4 * i + 2] = ilist->iatoms[2 * i + 1] + 1;
2136         ilist->iatoms[4 * i + 3] = ilist->iatoms[2 * i + 1] + 2;
2137     }
2138 }
2139
2140 static void do_ilists(gmx::ISerializer* serializer, InteractionLists* ilists, int file_version)
2141 {
2142     GMX_RELEASE_ASSERT(ilists, "Need a valid ilists object");
2143     GMX_RELEASE_ASSERT(ilists->size() == F_NRE,
2144                        "The code needs to be in sync with InteractionLists");
2145
2146     for (int j = 0; j < F_NRE; j++)
2147     {
2148         InteractionList& ilist  = (*ilists)[j];
2149         gmx_bool         bClear = FALSE;
2150         if (serializer->reading())
2151         {
2152             for (int k = 0; k < NFTUPD; k++)
2153             {
2154                 if ((file_version < ftupd[k].fvnr) && (j == ftupd[k].ftype))
2155                 {
2156                     bClear = TRUE;
2157                 }
2158             }
2159         }
2160         if (bClear)
2161         {
2162             ilist.iatoms.clear();
2163         }
2164         else
2165         {
2166             do_ilist(serializer, &ilist);
2167             if (file_version < 78 && j == F_SETTLE && !ilist.empty())
2168             {
2169                 add_settle_atoms(&ilist);
2170             }
2171         }
2172     }
2173 }
2174
2175 static void do_block(gmx::ISerializer* serializer, t_block* block)
2176 {
2177     serializer->doInt(&block->nr);
2178     if (serializer->reading())
2179     {
2180         if ((block->nalloc_index > 0) && (nullptr != block->index))
2181         {
2182             sfree(block->index);
2183         }
2184         block->nalloc_index = block->nr + 1;
2185         snew(block->index, block->nalloc_index);
2186     }
2187     serializer->doIntArray(block->index, block->nr + 1);
2188 }
2189
2190 static void doListOfLists(gmx::ISerializer* serializer, gmx::ListOfLists<int>* listOfLists)
2191 {
2192     int numLists = listOfLists->ssize();
2193     serializer->doInt(&numLists);
2194     int numElements = listOfLists->elementsView().ssize();
2195     serializer->doInt(&numElements);
2196     if (serializer->reading())
2197     {
2198         std::vector<int> listRanges(numLists + 1);
2199         serializer->doIntArray(listRanges.data(), numLists + 1);
2200         std::vector<int> elements(numElements);
2201         serializer->doIntArray(elements.data(), numElements);
2202         *listOfLists = gmx::ListOfLists<int>(std::move(listRanges), std::move(elements));
2203     }
2204     else
2205     {
2206         serializer->doIntArray(const_cast<int*>(listOfLists->listRangesView().data()), numLists + 1);
2207         serializer->doIntArray(const_cast<int*>(listOfLists->elementsView().data()), numElements);
2208     }
2209 }
2210
2211 /* This is a primitive routine to make it possible to translate atomic numbers
2212  * to element names when reading TPR files, without making the Gromacs library
2213  * directory a dependency on mdrun (which is the case if we need elements.dat).
2214  */
2215 static const char* atomicnumber_to_element(int atomicnumber)
2216 {
2217     const char* p;
2218
2219     /* This does not have to be complete, so we only include elements likely
2220      * to occur in PDB files.
2221      */
2222     switch (atomicnumber)
2223     {
2224         case 1: p = "H"; break;
2225         case 5: p = "B"; break;
2226         case 6: p = "C"; break;
2227         case 7: p = "N"; break;
2228         case 8: p = "O"; break;
2229         case 9: p = "F"; break;
2230         case 11: p = "Na"; break;
2231         case 12: p = "Mg"; break;
2232         case 15: p = "P"; break;
2233         case 16: p = "S"; break;
2234         case 17: p = "Cl"; break;
2235         case 18: p = "Ar"; break;
2236         case 19: p = "K"; break;
2237         case 20: p = "Ca"; break;
2238         case 25: p = "Mn"; break;
2239         case 26: p = "Fe"; break;
2240         case 28: p = "Ni"; break;
2241         case 29: p = "Cu"; break;
2242         case 30: p = "Zn"; break;
2243         case 35: p = "Br"; break;
2244         case 47: p = "Ag"; break;
2245         default: p = ""; break;
2246     }
2247     return p;
2248 }
2249
2250
2251 static void do_atom(gmx::ISerializer* serializer, t_atom* atom)
2252 {
2253     serializer->doReal(&atom->m);
2254     serializer->doReal(&atom->q);
2255     serializer->doReal(&atom->mB);
2256     serializer->doReal(&atom->qB);
2257     serializer->doUShort(&atom->type);
2258     serializer->doUShort(&atom->typeB);
2259     serializer->doEnumAsInt(&atom->ptype);
2260     serializer->doInt(&atom->resind);
2261     serializer->doInt(&atom->atomnumber);
2262     if (serializer->reading())
2263     {
2264         /* Set element string from atomic number if present.
2265          * This routine returns an empty string if the name is not found.
2266          */
2267         std::strncpy(atom->elem, atomicnumber_to_element(atom->atomnumber), 4);
2268         /* avoid warnings about potentially unterminated string */
2269         atom->elem[3] = '\0';
2270     }
2271 }
2272
2273 static void do_grps(gmx::ISerializer* serializer, gmx::ArrayRef<AtomGroupIndices> grps)
2274 {
2275     for (auto& group : grps)
2276     {
2277         int size = group.size();
2278         serializer->doInt(&size);
2279         if (serializer->reading())
2280         {
2281             group.resize(size);
2282         }
2283         serializer->doIntArray(group.data(), size);
2284     }
2285 }
2286
2287 static void do_symstr(gmx::ISerializer* serializer, char*** nm, t_symtab* symtab)
2288 {
2289     int ls;
2290
2291     if (serializer->reading())
2292     {
2293         serializer->doInt(&ls);
2294         *nm = get_symtab_handle(symtab, ls);
2295     }
2296     else
2297     {
2298         ls = lookup_symtab(symtab, *nm);
2299         serializer->doInt(&ls);
2300     }
2301 }
2302
2303 static void do_strstr(gmx::ISerializer* serializer, int nstr, char*** nm, t_symtab* symtab)
2304 {
2305     int j;
2306
2307     for (j = 0; (j < nstr); j++)
2308     {
2309         do_symstr(serializer, &(nm[j]), symtab);
2310     }
2311 }
2312
2313 static void do_resinfo(gmx::ISerializer* serializer, int n, t_resinfo* ri, t_symtab* symtab, int file_version)
2314 {
2315     int j;
2316
2317     for (j = 0; (j < n); j++)
2318     {
2319         do_symstr(serializer, &(ri[j].name), symtab);
2320         if (file_version >= 63)
2321         {
2322             serializer->doInt(&ri[j].nr);
2323             serializer->doUChar(&ri[j].ic);
2324         }
2325         else
2326         {
2327             ri[j].nr = j + 1;
2328             ri[j].ic = ' ';
2329         }
2330     }
2331 }
2332
2333 static void do_atoms(gmx::ISerializer* serializer, t_atoms* atoms, t_symtab* symtab, int file_version)
2334 {
2335     int i;
2336
2337     serializer->doInt(&atoms->nr);
2338     serializer->doInt(&atoms->nres);
2339     if (serializer->reading())
2340     {
2341         /* Since we have always written all t_atom properties in the tpr file
2342          * (at least for all backward compatible versions), we don't store
2343          * but simple set the booleans here.
2344          */
2345         atoms->haveMass    = TRUE;
2346         atoms->haveCharge  = TRUE;
2347         atoms->haveType    = TRUE;
2348         atoms->haveBState  = TRUE;
2349         atoms->havePdbInfo = FALSE;
2350
2351         snew(atoms->atom, atoms->nr);
2352         snew(atoms->atomname, atoms->nr);
2353         snew(atoms->atomtype, atoms->nr);
2354         snew(atoms->atomtypeB, atoms->nr);
2355         snew(atoms->resinfo, atoms->nres);
2356         atoms->pdbinfo = nullptr;
2357     }
2358     else
2359     {
2360         GMX_RELEASE_ASSERT(atoms->haveMass && atoms->haveCharge && atoms->haveType && atoms->haveBState,
2361                            "Mass, charge, atomtype and B-state parameters should be present in "
2362                            "t_atoms when writing a tpr file");
2363     }
2364     for (i = 0; (i < atoms->nr); i++)
2365     {
2366         do_atom(serializer, &atoms->atom[i]);
2367     }
2368     do_strstr(serializer, atoms->nr, atoms->atomname, symtab);
2369     do_strstr(serializer, atoms->nr, atoms->atomtype, symtab);
2370     do_strstr(serializer, atoms->nr, atoms->atomtypeB, symtab);
2371
2372     do_resinfo(serializer, atoms->nres, atoms->resinfo, symtab, file_version);
2373 }
2374
2375 static void do_groups(gmx::ISerializer* serializer, SimulationGroups* groups, t_symtab* symtab)
2376 {
2377     do_grps(serializer, groups->groups);
2378     int numberOfGroupNames = groups->groupNames.size();
2379     serializer->doInt(&numberOfGroupNames);
2380     if (serializer->reading())
2381     {
2382         groups->groupNames.resize(numberOfGroupNames);
2383     }
2384     do_strstr(serializer, numberOfGroupNames, groups->groupNames.data(), symtab);
2385     for (auto group : gmx::keysOf(groups->groupNumbers))
2386     {
2387         int numberOfGroupNumbers = groups->numberOfGroupNumbers(group);
2388         serializer->doInt(&numberOfGroupNumbers);
2389         if (numberOfGroupNumbers != 0)
2390         {
2391             if (serializer->reading())
2392             {
2393                 groups->groupNumbers[group].resize(numberOfGroupNumbers);
2394             }
2395             serializer->doUCharArray(groups->groupNumbers[group].data(), numberOfGroupNumbers);
2396         }
2397     }
2398 }
2399
2400 static void do_atomtypes(gmx::ISerializer* serializer, t_atomtypes* atomtypes, int file_version)
2401 {
2402     int j;
2403
2404     serializer->doInt(&atomtypes->nr);
2405     j = atomtypes->nr;
2406     if (serializer->reading())
2407     {
2408         snew(atomtypes->atomnumber, j);
2409     }
2410     if (serializer->reading() && file_version < tpxv_RemoveImplicitSolvation)
2411     {
2412         std::vector<real> dummy(atomtypes->nr, 0);
2413         serializer->doRealArray(dummy.data(), dummy.size());
2414         serializer->doRealArray(dummy.data(), dummy.size());
2415         serializer->doRealArray(dummy.data(), dummy.size());
2416     }
2417     serializer->doIntArray(atomtypes->atomnumber, j);
2418
2419     if (serializer->reading() && file_version >= 60 && file_version < tpxv_RemoveImplicitSolvation)
2420     {
2421         std::vector<real> dummy(atomtypes->nr, 0);
2422         serializer->doRealArray(dummy.data(), dummy.size());
2423         serializer->doRealArray(dummy.data(), dummy.size());
2424     }
2425 }
2426
2427 static void do_symtab(gmx::ISerializer* serializer, t_symtab* symtab)
2428 {
2429     int       i, nr;
2430     t_symbuf* symbuf;
2431
2432     serializer->doInt(&symtab->nr);
2433     nr = symtab->nr;
2434     if (serializer->reading())
2435     {
2436         snew(symtab->symbuf, 1);
2437         symbuf          = symtab->symbuf;
2438         symbuf->bufsize = nr;
2439         snew(symbuf->buf, nr);
2440         for (i = 0; (i < nr); i++)
2441         {
2442             std::string buf;
2443             serializer->doString(&buf);
2444             symbuf->buf[i] = gmx_strdup(buf.c_str());
2445         }
2446     }
2447     else
2448     {
2449         symbuf = symtab->symbuf;
2450         while (symbuf != nullptr)
2451         {
2452             for (i = 0; (i < symbuf->bufsize) && (i < nr); i++)
2453             {
2454                 std::string buf = symbuf->buf[i];
2455                 serializer->doString(&buf);
2456             }
2457             nr -= i;
2458             symbuf = symbuf->next;
2459         }
2460         if (nr != 0)
2461         {
2462             gmx_fatal(FARGS, "nr of symtab strings left: %d", nr);
2463         }
2464     }
2465 }
2466
2467 static void do_cmap(gmx::ISerializer* serializer, gmx_cmap_t* cmap_grid)
2468 {
2469
2470     int ngrid = cmap_grid->cmapdata.size();
2471     serializer->doInt(&ngrid);
2472     serializer->doInt(&cmap_grid->grid_spacing);
2473
2474     int gs    = cmap_grid->grid_spacing;
2475     int nelem = gs * gs;
2476
2477     if (serializer->reading())
2478     {
2479         cmap_grid->cmapdata.resize(ngrid);
2480
2481         for (int i = 0; i < ngrid; i++)
2482         {
2483             cmap_grid->cmapdata[i].cmap.resize(4 * nelem);
2484         }
2485     }
2486
2487     for (int i = 0; i < ngrid; i++)
2488     {
2489         for (int j = 0; j < nelem; j++)
2490         {
2491             serializer->doReal(&cmap_grid->cmapdata[i].cmap[j * 4]);
2492             serializer->doReal(&cmap_grid->cmapdata[i].cmap[j * 4 + 1]);
2493             serializer->doReal(&cmap_grid->cmapdata[i].cmap[j * 4 + 2]);
2494             serializer->doReal(&cmap_grid->cmapdata[i].cmap[j * 4 + 3]);
2495         }
2496     }
2497 }
2498
2499
2500 static void do_moltype(gmx::ISerializer* serializer, gmx_moltype_t* molt, t_symtab* symtab, int file_version)
2501 {
2502     do_symstr(serializer, &(molt->name), symtab);
2503
2504     do_atoms(serializer, &molt->atoms, symtab, file_version);
2505
2506     do_ilists(serializer, &molt->ilist, file_version);
2507
2508     /* TODO: Remove the obsolete charge group index from the file */
2509     t_block cgs;
2510     cgs.nr           = molt->atoms.nr;
2511     cgs.nalloc_index = cgs.nr + 1;
2512     snew(cgs.index, cgs.nalloc_index);
2513     for (int i = 0; i < cgs.nr + 1; i++)
2514     {
2515         cgs.index[i] = i;
2516     }
2517     do_block(serializer, &cgs);
2518     sfree(cgs.index);
2519
2520     /* This used to be in the atoms struct */
2521     doListOfLists(serializer, &molt->excls);
2522 }
2523
2524 static void do_molblock(gmx::ISerializer* serializer, gmx_molblock_t* molb, int numAtomsPerMolecule)
2525 {
2526     serializer->doInt(&molb->type);
2527     serializer->doInt(&molb->nmol);
2528     /* To maintain forward topology reading compatibility, we store #atoms.
2529      * TODO: Change this to conditional reading of a dummy int when we
2530      *       increase tpx_generation.
2531      */
2532     serializer->doInt(&numAtomsPerMolecule);
2533     /* Position restraint coordinates */
2534     int numPosres_xA = molb->posres_xA.size();
2535     serializer->doInt(&numPosres_xA);
2536     if (numPosres_xA > 0)
2537     {
2538         if (serializer->reading())
2539         {
2540             molb->posres_xA.resize(numPosres_xA);
2541         }
2542         serializer->doRvecArray(as_rvec_array(molb->posres_xA.data()), numPosres_xA);
2543     }
2544     int numPosres_xB = molb->posres_xB.size();
2545     serializer->doInt(&numPosres_xB);
2546     if (numPosres_xB > 0)
2547     {
2548         if (serializer->reading())
2549         {
2550             molb->posres_xB.resize(numPosres_xB);
2551         }
2552         serializer->doRvecArray(as_rvec_array(molb->posres_xB.data()), numPosres_xB);
2553     }
2554 }
2555
2556 static void set_disres_npair(gmx_mtop_t* mtop)
2557 {
2558     gmx::ArrayRef<t_iparams> ip = mtop->ffparams.iparams;
2559
2560     for (const auto ilist : IListRange(*mtop))
2561     {
2562         const InteractionList& il = ilist.list()[F_DISRES];
2563
2564         if (!il.empty())
2565         {
2566             gmx::ArrayRef<const int> a     = il.iatoms;
2567             int                      npair = 0;
2568             for (int i = 0; i < il.size(); i += 3)
2569             {
2570                 npair++;
2571                 if (i + 3 == il.size() || ip[a[i]].disres.label != ip[a[i + 3]].disres.label)
2572                 {
2573                     ip[a[i]].disres.npair = npair;
2574                     npair                 = 0;
2575                 }
2576             }
2577         }
2578     }
2579 }
2580
2581 static void do_mtop(gmx::ISerializer* serializer, gmx_mtop_t* mtop, int file_version)
2582 {
2583     do_symtab(serializer, &(mtop->symtab));
2584
2585     do_symstr(serializer, &(mtop->name), &(mtop->symtab));
2586
2587     do_ffparams(serializer, &mtop->ffparams, file_version);
2588
2589     int nmoltype = mtop->moltype.size();
2590     serializer->doInt(&nmoltype);
2591     if (serializer->reading())
2592     {
2593         mtop->moltype.resize(nmoltype);
2594     }
2595     for (gmx_moltype_t& moltype : mtop->moltype)
2596     {
2597         do_moltype(serializer, &moltype, &mtop->symtab, file_version);
2598     }
2599
2600     int nmolblock = mtop->molblock.size();
2601     serializer->doInt(&nmolblock);
2602     if (serializer->reading())
2603     {
2604         mtop->molblock.resize(nmolblock);
2605     }
2606     for (gmx_molblock_t& molblock : mtop->molblock)
2607     {
2608         int numAtomsPerMolecule = (serializer->reading() ? 0 : mtop->moltype[molblock.type].atoms.nr);
2609         do_molblock(serializer, &molblock, numAtomsPerMolecule);
2610     }
2611     serializer->doInt(&mtop->natoms);
2612
2613     if (file_version >= tpxv_IntermolecularBondeds)
2614     {
2615         serializer->doBool(&mtop->bIntermolecularInteractions);
2616         if (mtop->bIntermolecularInteractions)
2617         {
2618             if (serializer->reading())
2619             {
2620                 mtop->intermolecular_ilist = std::make_unique<InteractionLists>();
2621             }
2622             do_ilists(serializer, mtop->intermolecular_ilist.get(), file_version);
2623         }
2624     }
2625     else
2626     {
2627         mtop->bIntermolecularInteractions = FALSE;
2628     }
2629
2630     do_atomtypes(serializer, &(mtop->atomtypes), file_version);
2631
2632     if (file_version >= 65)
2633     {
2634         do_cmap(serializer, &mtop->ffparams.cmap_grid);
2635     }
2636     else
2637     {
2638         mtop->ffparams.cmap_grid.grid_spacing = 0;
2639         mtop->ffparams.cmap_grid.cmapdata.clear();
2640     }
2641
2642     do_groups(serializer, &mtop->groups, &(mtop->symtab));
2643
2644     mtop->haveMoleculeIndices = true;
2645
2646     if (file_version >= tpxv_StoreNonBondedInteractionExclusionGroup)
2647     {
2648         std::int64_t intermolecularExclusionGroupSize = gmx::ssize(mtop->intermolecularExclusionGroup);
2649         serializer->doInt64(&intermolecularExclusionGroupSize);
2650         GMX_RELEASE_ASSERT(intermolecularExclusionGroupSize >= 0,
2651                            "Number of atoms with excluded intermolecular non-bonded interactions "
2652                            "is negative.");
2653         mtop->intermolecularExclusionGroup.resize(intermolecularExclusionGroupSize); // no effect when writing
2654         serializer->doIntArray(mtop->intermolecularExclusionGroup.data(),
2655                                mtop->intermolecularExclusionGroup.size());
2656     }
2657
2658     if (serializer->reading())
2659     {
2660         close_symtab(&(mtop->symtab));
2661     }
2662 }
2663
2664 /*! \brief
2665  * Read the first part of the TPR file to find general system information.
2666  *
2667  * If \p TopOnlyOK is true then we can read even future versions
2668  * of tpx files, provided the \p fileGeneration hasn't changed.
2669  * If it is false, we need the \p ir too, and bail out
2670  * if the file is newer than the program.
2671  *
2672  * The version and generation of the topology (see top of this file)
2673  * are returned in the two last arguments, if those arguments are non-nullptr.
2674  *
2675  * If possible, we will read the \p ir even when \p TopOnlyOK is true.
2676  *
2677  * \param[in,out] serializer The serializer used to handle header processing.
2678  * \param[in,out] tpx File header datastructure.
2679  * \param[in]     filename The name of the file being read/written
2680  * \param[in,out] fio File handle.
2681  * \param[in] TopOnlyOK If not reading \p ir is fine or not.
2682  */
2683 static void do_tpxheader(gmx::FileIOXdrSerializer* serializer,
2684                          TpxFileHeader*            tpx,
2685                          const char*               filename,
2686                          t_fileio*                 fio,
2687                          bool                      TopOnlyOK)
2688 {
2689     int  precision;
2690     int  idum = 0;
2691     real rdum = 0;
2692
2693     /* XDR binary topology file */
2694     precision = sizeof(real);
2695     std::string buf;
2696     std::string fileTag;
2697     if (serializer->reading())
2698     {
2699         serializer->doString(&buf);
2700         if (std::strncmp(buf.c_str(), "VERSION", 7) != 0)
2701         {
2702             gmx_fatal(
2703                     FARGS,
2704                     "Can not read file %s,\n"
2705                     "             this file is from a GROMACS version which is older than 2.0\n"
2706                     "             Make a new one with grompp or use a gro or pdb file, if possible",
2707                     filename);
2708         }
2709         // We need to know the precision used to write the TPR file, to match it
2710         // to the precision of the currently running binary. If the precisions match
2711         // there is no problem, but mismatching precision needs to be accounted for
2712         // by reading into temporary variables of the correct precision instead
2713         // of the desired target datastructures.
2714         serializer->doInt(&precision);
2715         tpx->isDouble = (precision == sizeof(double));
2716         if ((precision != sizeof(float)) && !tpx->isDouble)
2717         {
2718             gmx_fatal(FARGS,
2719                       "Unknown precision in file %s: real is %d bytes "
2720                       "instead of %zu or %zu",
2721                       filename,
2722                       precision,
2723                       sizeof(float),
2724                       sizeof(double));
2725         }
2726         gmx_fio_setprecision(fio, tpx->isDouble);
2727         fprintf(stderr,
2728                 "Reading file %s, %s (%s precision)\n",
2729                 filename,
2730                 buf.c_str(),
2731                 tpx->isDouble ? "double" : "single");
2732     }
2733     else
2734     {
2735         buf = gmx::formatString("VERSION %s", gmx_version());
2736         serializer->doString(&buf);
2737         gmx_fio_setprecision(fio, tpx->isDouble);
2738         serializer->doInt(&precision);
2739         fileTag = tpx_tag;
2740     }
2741
2742     /* Check versions! */
2743     serializer->doInt(&tpx->fileVersion);
2744
2745     /* This is for backward compatibility with development versions 77-79
2746      * where the tag was, mistakenly, placed before the generation,
2747      * which would cause a segv instead of a proper error message
2748      * when reading the topology only from tpx with <77 code.
2749      */
2750     if (tpx->fileVersion >= 77 && tpx->fileVersion <= 79)
2751     {
2752         serializer->doString(&fileTag);
2753     }
2754
2755     serializer->doInt(&tpx->fileGeneration);
2756
2757     if (tpx->fileVersion >= 81)
2758     {
2759         serializer->doString(&fileTag);
2760     }
2761     if (serializer->reading())
2762     {
2763         if (tpx->fileVersion < 77)
2764         {
2765             /* Versions before 77 don't have the tag, set it to release */
2766             fileTag = gmx::formatString("%s", TPX_TAG_RELEASE);
2767         }
2768
2769         if (fileTag != tpx_tag)
2770         {
2771             fprintf(stderr, "Note: file tpx tag '%s', software tpx tag '%s'\n", fileTag.c_str(), tpx_tag.c_str());
2772
2773             /* We only support reading tpx files with the same tag as the code
2774              * or tpx files with the release tag and with lower version number.
2775              */
2776             if (fileTag != TPX_TAG_RELEASE && tpx->fileVersion < tpx_version)
2777             {
2778                 gmx_fatal(FARGS,
2779                           "tpx tag/version mismatch: reading tpx file (%s) version %d, tag '%s' "
2780                           "with program for tpx version %d, tag '%s'",
2781                           filename,
2782                           tpx->fileVersion,
2783                           fileTag.c_str(),
2784                           tpx_version,
2785                           tpx_tag.c_str());
2786             }
2787         }
2788     }
2789
2790     if ((tpx->fileVersion <= tpx_incompatible_version)
2791         || ((tpx->fileVersion > tpx_version) && !TopOnlyOK) || (tpx->fileGeneration > tpx_generation)
2792         || tpx_version == 80) /*80 was used by both 5.0-dev and 4.6-dev*/
2793     {
2794         gmx_fatal(FARGS,
2795                   "reading tpx file (%s) version %d with version %d program",
2796                   filename,
2797                   tpx->fileVersion,
2798                   tpx_version);
2799     }
2800
2801     serializer->doInt(&tpx->natoms);
2802     serializer->doInt(&tpx->ngtc);
2803
2804     if (tpx->fileVersion < 62)
2805     {
2806         serializer->doInt(&idum);
2807         serializer->doReal(&rdum);
2808     }
2809     if (tpx->fileVersion >= 79)
2810     {
2811         serializer->doInt(&tpx->fep_state);
2812     }
2813     serializer->doReal(&tpx->lambda);
2814     serializer->doBool(&tpx->bIr);
2815     serializer->doBool(&tpx->bTop);
2816     serializer->doBool(&tpx->bX);
2817     serializer->doBool(&tpx->bV);
2818     serializer->doBool(&tpx->bF);
2819     serializer->doBool(&tpx->bBox);
2820
2821     if (tpx->fileVersion >= tpxv_AddSizeField && tpx->fileGeneration >= 27)
2822     {
2823         if (!serializer->reading())
2824         {
2825             GMX_RELEASE_ASSERT(tpx->sizeOfTprBody != 0,
2826                                "Not possible to write new file with zero TPR body size");
2827         }
2828         serializer->doInt64(&tpx->sizeOfTprBody);
2829     }
2830
2831     if ((tpx->fileGeneration > tpx_generation))
2832     {
2833         /* This can only happen if TopOnlyOK=TRUE */
2834         tpx->bIr = FALSE;
2835     }
2836 }
2837
2838 #define do_test(serializer, b, p)                            \
2839     if ((serializer)->reading() && ((p) != nullptr) && !(b)) \
2840     gmx_fatal(FARGS, "No %s in input file", #p)
2841
2842 /*! \brief
2843  * Process the first part of the TPR into the state datastructure.
2844  *
2845  * Due to the structure of the legacy code, it is necessary
2846  * to split up the state reading into two parts, with the
2847  * box and legacy temperature coupling processed before the
2848  * topology datastructures.
2849  *
2850  * See the documentation for do_tpx_body for the correct order of
2851  * the operations for reading a tpr file.
2852  *
2853  * \param[in] serializer Abstract serializer used to read/write data.
2854  * \param[in] tpx The file header data.
2855  * \param[in, out] state Global state data.
2856  */
2857 static void do_tpx_state_first(gmx::ISerializer* serializer, TpxFileHeader* tpx, t_state* state)
2858 {
2859     if (serializer->reading())
2860     {
2861         state->flags = 0;
2862         init_gtc_state(state, tpx->ngtc, 0, 0);
2863     }
2864     do_test(serializer, tpx->bBox, state->box);
2865     if (tpx->bBox)
2866     {
2867         serializer->doRvecArray(state->box, DIM);
2868         if (tpx->fileVersion >= 51)
2869         {
2870             serializer->doRvecArray(state->box_rel, DIM);
2871         }
2872         else
2873         {
2874             /* We initialize box_rel after reading the inputrec */
2875             clear_mat(state->box_rel);
2876         }
2877         serializer->doRvecArray(state->boxv, DIM);
2878         if (tpx->fileVersion < 56)
2879         {
2880             matrix mdum;
2881             serializer->doRvecArray(mdum, DIM);
2882         }
2883     }
2884
2885     if (state->ngtc > 0)
2886     {
2887         real* dumv;
2888         snew(dumv, state->ngtc);
2889         if (tpx->fileVersion < 69)
2890         {
2891             serializer->doRealArray(dumv, state->ngtc);
2892         }
2893         /* These used to be the Berendsen tcoupl_lambda's */
2894         serializer->doRealArray(dumv, state->ngtc);
2895         sfree(dumv);
2896     }
2897 }
2898
2899 /*! \brief
2900  * Process global topology data.
2901  *
2902  * See the documentation for do_tpx_body for the correct order of
2903  * the operations for reading a tpr file.
2904  *
2905  * \param[in] serializer Abstract serializer  used to read/write data.
2906  * \param[in] tpx The file header data.
2907  * \param[in,out] mtop Global topology.
2908  */
2909 static void do_tpx_mtop(gmx::ISerializer* serializer, TpxFileHeader* tpx, gmx_mtop_t* mtop)
2910 {
2911     do_test(serializer, tpx->bTop, mtop);
2912     if (tpx->bTop)
2913     {
2914         if (mtop)
2915         {
2916             do_mtop(serializer, mtop, tpx->fileVersion);
2917             set_disres_npair(mtop);
2918             mtop->finalize();
2919         }
2920         else
2921         {
2922             gmx_mtop_t dum_top;
2923             do_mtop(serializer, &dum_top, tpx->fileVersion);
2924         }
2925     }
2926 }
2927 /*! \brief
2928  * Process coordinate vectors for state data.
2929  *
2930  * Main part of state gets processed here.
2931  *
2932  * See the documentation for do_tpx_body for the correct order of
2933  * the operations for reading a tpr file.
2934  *
2935  * \param[in] serializer Abstract serializer used to read/write data.
2936  * \param[in] tpx The file header data.
2937  * \param[in,out] state Global state data.
2938  * \param[in,out] x Individual coordinates for processing, deprecated.
2939  * \param[in,out] v Individual velocities for processing, deprecated.
2940  */
2941 static void do_tpx_state_second(gmx::ISerializer* serializer, TpxFileHeader* tpx, t_state* state, rvec* x, rvec* v)
2942 {
2943     if (!serializer->reading())
2944     {
2945         GMX_RELEASE_ASSERT(
2946                 x == nullptr && v == nullptr,
2947                 "Passing separate x and v pointers to do_tpx() is not supported when writing");
2948     }
2949     else
2950     {
2951         GMX_RELEASE_ASSERT(!(x == nullptr && v != nullptr),
2952                            "Passing x==NULL and v!=NULL is not supported");
2953     }
2954
2955     if (serializer->reading())
2956     {
2957         if (x == nullptr)
2958         {
2959             // v is also nullptr by the above assertion, so we may
2960             // need to make memory in state for storing the contents
2961             // of the tpx file.
2962             if (tpx->bX)
2963             {
2964                 state->flags |= enumValueToBitMask(StateEntry::X);
2965             }
2966             if (tpx->bV)
2967             {
2968                 state->flags |= enumValueToBitMask(StateEntry::V);
2969             }
2970             state_change_natoms(state, tpx->natoms);
2971         }
2972     }
2973
2974     if (x == nullptr)
2975     {
2976         x = state->x.rvec_array();
2977         v = state->v.rvec_array();
2978     }
2979     do_test(serializer, tpx->bX, x);
2980     if (tpx->bX)
2981     {
2982         if (serializer->reading())
2983         {
2984             state->flags |= enumValueToBitMask(StateEntry::X);
2985         }
2986         serializer->doRvecArray(x, tpx->natoms);
2987     }
2988
2989     do_test(serializer, tpx->bV, v);
2990     if (tpx->bV)
2991     {
2992         if (serializer->reading())
2993         {
2994             state->flags |= enumValueToBitMask(StateEntry::V);
2995         }
2996         if (!v)
2997         {
2998             std::vector<gmx::RVec> dummyVelocities(tpx->natoms);
2999             serializer->doRvecArray(as_rvec_array(dummyVelocities.data()), tpx->natoms);
3000         }
3001         else
3002         {
3003             serializer->doRvecArray(v, tpx->natoms);
3004         }
3005     }
3006
3007     // No need to run do_test when the last argument is NULL
3008     if (tpx->bF)
3009     {
3010         std::vector<gmx::RVec> dummyForces(state->natoms);
3011         serializer->doRvecArray(as_rvec_array(dummyForces.data()), tpx->natoms);
3012     }
3013 }
3014 /*! \brief
3015  * Process simulation parameters.
3016  *
3017  * See the documentation for do_tpx_body for the correct order of
3018  * the operations for reading a tpr file.
3019  *
3020  * \param[in] serializer Abstract serializer used to read/write data.
3021  * \param[in] tpx The file header data.
3022  * \param[in,out] ir Datastructure with simulation parameters.
3023  */
3024 static PbcType do_tpx_ir(gmx::ISerializer* serializer, TpxFileHeader* tpx, t_inputrec* ir)
3025 {
3026     PbcType  pbcType;
3027     gmx_bool bPeriodicMols;
3028
3029     /* Starting with tpx version 26, we have the inputrec
3030      * at the end of the file, so we can ignore it
3031      * if the file is never than the software (but still the
3032      * same generation - see comments at the top of this file.
3033      *
3034      *
3035      */
3036     pbcType       = PbcType::Unset;
3037     bPeriodicMols = FALSE;
3038
3039     do_test(serializer, tpx->bIr, ir);
3040     if (tpx->bIr)
3041     {
3042         if (tpx->fileVersion >= 53)
3043         {
3044             /* Removed the pbc info from do_inputrec, since we always want it */
3045             if (!serializer->reading())
3046             {
3047                 pbcType       = ir->pbcType;
3048                 bPeriodicMols = ir->bPeriodicMols;
3049             }
3050             serializer->doInt(reinterpret_cast<int*>(&pbcType));
3051             serializer->doBool(&bPeriodicMols);
3052         }
3053         if (tpx->fileGeneration <= tpx_generation && ir)
3054         {
3055             do_inputrec(serializer, ir, tpx->fileVersion);
3056             if (tpx->fileVersion < 53)
3057             {
3058                 pbcType       = ir->pbcType;
3059                 bPeriodicMols = ir->bPeriodicMols;
3060             }
3061         }
3062         if (serializer->reading() && ir && tpx->fileVersion >= 53)
3063         {
3064             /* We need to do this after do_inputrec, since that initializes ir */
3065             ir->pbcType       = pbcType;
3066             ir->bPeriodicMols = bPeriodicMols;
3067         }
3068     }
3069     return pbcType;
3070 }
3071
3072 /*! \brief
3073  * Correct and finalize read information.
3074  *
3075  * If \p state is nullptr, skip the parts dependent on it.
3076  *
3077  * See the documentation for do_tpx_body for the correct order of
3078  * the operations for reading a tpr file.
3079  *
3080  * \param[in] tpx The file header used to check version numbers.
3081  * \param[out] ir Input rec that needs correction.
3082  * \param[out] state State needing correction.
3083  * \param[out] mtop Topology to finalize.
3084  */
3085 static void do_tpx_finalize(TpxFileHeader* tpx, t_inputrec* ir, t_state* state, gmx_mtop_t* mtop)
3086 {
3087     if (tpx->fileVersion < 51 && state)
3088     {
3089         set_box_rel(ir, state);
3090     }
3091     if (tpx->bIr && ir)
3092     {
3093         if (state && state->ngtc == 0)
3094         {
3095             /* Reading old version without tcoupl state data: set it */
3096             init_gtc_state(state, ir->opts.ngtc, 0, ir->opts.nhchainlength);
3097         }
3098         if (tpx->bTop && mtop)
3099         {
3100             if (tpx->fileVersion < 57)
3101             {
3102                 ir->eDisre = !mtop->moltype[0].ilist[F_DISRES].empty()
3103                                      ? DistanceRestraintRefinement::Simple
3104                                      : DistanceRestraintRefinement::None;
3105             }
3106         }
3107     }
3108 }
3109
3110 /*! \brief
3111  * Process TPR data for file reading/writing.
3112  *
3113  * The TPR file gets processed in in four stages due to the organization
3114  * of the data within it.
3115  *
3116  * First, state data for the box is processed in do_tpx_state_first.
3117  * This is followed by processing the topology in do_tpx_mtop.
3118  * Coordinate and velocity vectors are handled next in do_tpx_state_second.
3119  * The last file information processed is the collection of simulation parameters in do_tpx_ir.
3120  * When reading, a final processing step is undertaken at the end.
3121  *
3122  * \param[in] serializer Abstract serializer used to read/write data.
3123  * \param[in] tpx The file header data.
3124  * \param[in,out] ir Datastructures with simulation parameters.
3125  * \param[in,out] state Global state data.
3126  * \param[in,out] x Individual coordinates for processing, deprecated.
3127  * \param[in,out] v Individual velocities for processing, deprecated.
3128  * \param[in,out] mtop Global topology.
3129  */
3130 static PbcType do_tpx_body(gmx::ISerializer* serializer,
3131                            TpxFileHeader*    tpx,
3132                            t_inputrec*       ir,
3133                            t_state*          state,
3134                            rvec*             x,
3135                            rvec*             v,
3136                            gmx_mtop_t*       mtop)
3137 {
3138     if (state)
3139     {
3140         do_tpx_state_first(serializer, tpx, state);
3141     }
3142     do_tpx_mtop(serializer, tpx, mtop);
3143     if (state)
3144     {
3145         do_tpx_state_second(serializer, tpx, state, x, v);
3146     }
3147     PbcType pbcType = do_tpx_ir(serializer, tpx, ir);
3148     if (serializer->reading())
3149     {
3150         do_tpx_finalize(tpx, ir, state, mtop);
3151     }
3152     return pbcType;
3153 }
3154
3155 /*! \brief
3156  * Overload for do_tpx_body that defaults to state vectors being nullptr.
3157  *
3158  * \param[in] serializer Abstract serializer used to read/write data.
3159  * \param[in] tpx The file header data.
3160  * \param[in,out] ir Datastructures with simulation parameters.
3161  * \param[in,out] mtop Global topology.
3162  */
3163 static PbcType do_tpx_body(gmx::ISerializer* serializer, TpxFileHeader* tpx, t_inputrec* ir, gmx_mtop_t* mtop)
3164 {
3165     return do_tpx_body(serializer, tpx, ir, nullptr, nullptr, nullptr, mtop);
3166 }
3167
3168 static t_fileio* open_tpx(const char* fn, const char* mode)
3169 {
3170     return gmx_fio_open(fn, mode);
3171 }
3172
3173 static void close_tpx(t_fileio* fio)
3174 {
3175     gmx_fio_close(fio);
3176 }
3177
3178 /*! \brief
3179  * Fill information into the header only from state before writing.
3180  *
3181  * Populating the header needs to be independent from writing the information
3182  * to file to allow things like writing the raw byte stream.
3183  *
3184  * \param[in] state The current simulation state. Can't write without it.
3185  * \param[in] ir Parameter and system information.
3186  * \param[in] mtop Global topology.
3187  * \returns Fully populated header.
3188  */
3189 static TpxFileHeader populateTpxHeader(const t_state& state, const t_inputrec* ir, const gmx_mtop_t* mtop)
3190 {
3191     TpxFileHeader header;
3192     header.natoms         = state.natoms;
3193     header.ngtc           = state.ngtc;
3194     header.fep_state      = state.fep_state;
3195     header.lambda         = state.lambda[FreeEnergyPerturbationCouplingType::Fep];
3196     header.bIr            = ir != nullptr;
3197     header.bTop           = mtop != nullptr;
3198     header.bX             = (state.flags & enumValueToBitMask(StateEntry::X)) != 0;
3199     header.bV             = (state.flags & enumValueToBitMask(StateEntry::V)) != 0;
3200     header.bF             = false;
3201     header.bBox           = true;
3202     header.fileVersion    = tpx_version;
3203     header.fileGeneration = tpx_generation;
3204     header.isDouble       = (sizeof(real) == sizeof(double));
3205
3206     return header;
3207 }
3208
3209 /*! \brief
3210  * Process the body of a TPR file as an opaque data buffer.
3211  *
3212  * Reads/writes the information in \p buffer from/to the \p serializer
3213  * provided to the function. Does not interact with the actual
3214  * TPR datastructures but with an in memory representation of the
3215  * data, so that this data can be efficiently read or written from/to
3216  * an original source.
3217  *
3218  * \param[in] serializer The abstract serializer used for reading or writing
3219  *                       the information in \p buffer.
3220  * \param[in,out] buffer Information from TPR file as char buffer.
3221  */
3222 static void doTpxBodyBuffer(gmx::ISerializer* serializer, gmx::ArrayRef<char> buffer)
3223 {
3224     serializer->doOpaque(buffer.data(), buffer.size());
3225 }
3226
3227 /*! \brief
3228  * Populates simulation datastructures.
3229  *
3230  * Here the information from the serialization interface \p serializer
3231  * is used to first populate the datastructures containing the simulation
3232  * information. Depending on the version found in the header \p tpx,
3233  * this is done using the new reading of the data as one block from disk,
3234  * followed by complete deserialization of the information read from there.
3235  * Otherwise, the datastructures are populated as before one by one from disk.
3236  * The second version is the default for the legacy tools that read the
3237  * coordinates and velocities separate from the state.
3238  *
3239  * After reading in the data, a separate buffer is populated from them
3240  * containing only \p ir and \p mtop that can be communicated directly
3241  * to nodes needing the information to set up a simulation.
3242  *
3243  * \param[in] tpx The file header.
3244  * \param[in] serializer The Serialization interface used to read the TPR.
3245  * \param[out] ir Input rec to populate.
3246  * \param[out] state State vectors to populate.
3247  * \param[out] x Coordinates to populate if needed.
3248  * \param[out] v Velocities to populate if needed.
3249  * \param[out] mtop Global topology to populate.
3250  *
3251  * \returns Partial de-serialized TPR used for communication to nodes.
3252  */
3253 static PartialDeserializedTprFile readTpxBody(TpxFileHeader*    tpx,
3254                                               gmx::ISerializer* serializer,
3255                                               t_inputrec*       ir,
3256                                               t_state*          state,
3257                                               rvec*             x,
3258                                               rvec*             v,
3259                                               gmx_mtop_t*       mtop)
3260 {
3261     PartialDeserializedTprFile partialDeserializedTpr;
3262     if (tpx->fileVersion >= tpxv_AddSizeField && tpx->fileGeneration >= 27)
3263     {
3264         partialDeserializedTpr.body.resize(tpx->sizeOfTprBody);
3265         partialDeserializedTpr.header = *tpx;
3266         doTpxBodyBuffer(serializer, partialDeserializedTpr.body);
3267
3268         partialDeserializedTpr.pbcType =
3269                 completeTprDeserialization(&partialDeserializedTpr, ir, state, x, v, mtop);
3270     }
3271     else
3272     {
3273         partialDeserializedTpr.pbcType = do_tpx_body(serializer, tpx, ir, state, x, v, mtop);
3274     }
3275     // Update header to system info for communication to nodes.
3276     // As we only need to communicate the inputrec and mtop to other nodes,
3277     // we prepare a new char buffer with the information we have already read
3278     // in on master.
3279     partialDeserializedTpr.header = populateTpxHeader(*state, ir, mtop);
3280     // Long-term we should move to use little endian in files to avoid extra byte swapping,
3281     // but since we just used the default XDR format (which is big endian) for the TPR
3282     // header it would cause third-party libraries reading our raw data to tear their hair
3283     // if we swap the endian in the middle of the file, so we stick to big endian in the
3284     // TPR file for now - and thus we ask the serializer to swap if this host is little endian.
3285     gmx::InMemorySerializer tprBodySerializer(gmx::EndianSwapBehavior::SwapIfHostIsLittleEndian);
3286     do_tpx_body(&tprBodySerializer, &partialDeserializedTpr.header, ir, mtop);
3287     partialDeserializedTpr.body = tprBodySerializer.finishAndGetBuffer();
3288
3289     return partialDeserializedTpr;
3290 }
3291
3292 /************************************************************
3293  *
3294  *  The following routines are the exported ones
3295  *
3296  ************************************************************/
3297
3298 TpxFileHeader readTpxHeader(const char* fileName, bool canReadTopologyOnly)
3299 {
3300     t_fileio* fio;
3301
3302     fio = open_tpx(fileName, "r");
3303     gmx::FileIOXdrSerializer serializer(fio);
3304
3305     TpxFileHeader tpx;
3306     do_tpxheader(&serializer, &tpx, fileName, fio, canReadTopologyOnly);
3307     close_tpx(fio);
3308     return tpx;
3309 }
3310
3311 void write_tpx_state(const char* fn, const t_inputrec* ir, const t_state* state, const gmx_mtop_t& mtop)
3312 {
3313     /* To write a state, we first need to write the state information to a buffer before
3314      * we append the raw bytes to the file. For this, the header information needs to be
3315      * populated before we write the main body because it has some information that is
3316      * otherwise not available.
3317      */
3318
3319     t_fileio* fio;
3320
3321     TpxFileHeader tpx = populateTpxHeader(*state, ir, &mtop);
3322     // Long-term we should move to use little endian in files to avoid extra byte swapping,
3323     // but since we just used the default XDR format (which is big endian) for the TPR
3324     // header it would cause third-party libraries reading our raw data to tear their hair
3325     // if we swap the endian in the middle of the file, so we stick to big endian in the
3326     // TPR file for now - and thus we ask the serializer to swap if this host is little endian.
3327     gmx::InMemorySerializer tprBodySerializer(gmx::EndianSwapBehavior::SwapIfHostIsLittleEndian);
3328
3329     do_tpx_body(&tprBodySerializer,
3330                 &tpx,
3331                 const_cast<t_inputrec*>(ir),
3332                 const_cast<t_state*>(state),
3333                 nullptr,
3334                 nullptr,
3335                 const_cast<gmx_mtop_t*>(&mtop));
3336
3337     std::vector<char> tprBody = tprBodySerializer.finishAndGetBuffer();
3338     tpx.sizeOfTprBody         = tprBody.size();
3339
3340     fio = open_tpx(fn, "w");
3341     gmx::FileIOXdrSerializer serializer(fio);
3342     do_tpxheader(&serializer, &tpx, fn, fio, ir == nullptr);
3343     doTpxBodyBuffer(&serializer, tprBody);
3344
3345     close_tpx(fio);
3346 }
3347
3348 PbcType completeTprDeserialization(PartialDeserializedTprFile* partialDeserializedTpr,
3349                                    t_inputrec*                 ir,
3350                                    t_state*                    state,
3351                                    rvec*                       x,
3352                                    rvec*                       v,
3353                                    gmx_mtop_t*                 mtop)
3354 {
3355     // Long-term we should move to use little endian in files to avoid extra byte swapping,
3356     // but since we just used the default XDR format (which is big endian) for the TPR
3357     // header it would cause third-party libraries reading our raw data to tear their hair
3358     // if we swap the endian in the middle of the file, so we stick to big endian in the
3359     // TPR file for now - and thus we ask the serializer to swap if this host is little endian.
3360     gmx::InMemoryDeserializer tprBodyDeserializer(partialDeserializedTpr->body,
3361                                                   partialDeserializedTpr->header.isDouble,
3362                                                   gmx::EndianSwapBehavior::SwapIfHostIsLittleEndian);
3363     return do_tpx_body(&tprBodyDeserializer, &partialDeserializedTpr->header, ir, state, x, v, mtop);
3364 }
3365
3366 PbcType completeTprDeserialization(PartialDeserializedTprFile* partialDeserializedTpr,
3367                                    t_inputrec*                 ir,
3368                                    gmx_mtop_t*                 mtop)
3369 {
3370     return completeTprDeserialization(partialDeserializedTpr, ir, nullptr, nullptr, nullptr, mtop);
3371 }
3372
3373 PartialDeserializedTprFile read_tpx_state(const char* fn, t_inputrec* ir, t_state* state, gmx_mtop_t* mtop)
3374 {
3375     t_fileio* fio;
3376     fio = open_tpx(fn, "r");
3377     gmx::FileIOXdrSerializer   serializer(fio);
3378     PartialDeserializedTprFile partialDeserializedTpr;
3379     do_tpxheader(&serializer, &partialDeserializedTpr.header, fn, fio, ir == nullptr);
3380     partialDeserializedTpr =
3381             readTpxBody(&partialDeserializedTpr.header, &serializer, ir, state, nullptr, nullptr, mtop);
3382     close_tpx(fio);
3383     return partialDeserializedTpr;
3384 }
3385
3386 PbcType read_tpx(const char* fn, t_inputrec* ir, matrix box, int* natoms, rvec* x, rvec* v, gmx_mtop_t* mtop)
3387 {
3388     t_fileio* fio;
3389     t_state   state;
3390
3391     TpxFileHeader tpx;
3392     fio = open_tpx(fn, "r");
3393     gmx::FileIOXdrSerializer serializer(fio);
3394     do_tpxheader(&serializer, &tpx, fn, fio, ir == nullptr);
3395     PartialDeserializedTprFile partialDeserializedTpr =
3396             readTpxBody(&tpx, &serializer, ir, &state, x, v, mtop);
3397     close_tpx(fio);
3398     if (mtop != nullptr && natoms != nullptr)
3399     {
3400         *natoms = mtop->natoms;
3401     }
3402     if (box)
3403     {
3404         copy_mat(state.box, box);
3405     }
3406     return partialDeserializedTpr.pbcType;
3407 }
3408
3409 PbcType read_tpx_top(const char* fn, t_inputrec* ir, matrix box, int* natoms, rvec* x, rvec* v, t_topology* top)
3410 {
3411     gmx_mtop_t mtop;
3412     PbcType    pbcType;
3413
3414     pbcType = read_tpx(fn, ir, box, natoms, x, v, &mtop);
3415
3416     *top = gmx_mtop_t_to_t_topology(&mtop, true);
3417
3418     return pbcType;
3419 }
3420
3421 gmx_bool fn2bTPX(const char* file)
3422 {
3423     return (efTPR == fn2ftp(file));
3424 }
3425
3426 void pr_tpxheader(FILE* fp, int indent, const char* title, const TpxFileHeader* sh)
3427 {
3428     if (available(fp, sh, indent, title))
3429     {
3430         indent = pr_title(fp, indent, title);
3431         pr_indent(fp, indent);
3432         fprintf(fp, "bIr    = %spresent\n", sh->bIr ? "" : "not ");
3433         pr_indent(fp, indent);
3434         fprintf(fp, "bBox   = %spresent\n", sh->bBox ? "" : "not ");
3435         pr_indent(fp, indent);
3436         fprintf(fp, "bTop   = %spresent\n", sh->bTop ? "" : "not ");
3437         pr_indent(fp, indent);
3438         fprintf(fp, "bX     = %spresent\n", sh->bX ? "" : "not ");
3439         pr_indent(fp, indent);
3440         fprintf(fp, "bV     = %spresent\n", sh->bV ? "" : "not ");
3441         pr_indent(fp, indent);
3442         fprintf(fp, "bF     = %spresent\n", sh->bF ? "" : "not ");
3443
3444         pr_indent(fp, indent);
3445         fprintf(fp, "natoms = %d\n", sh->natoms);
3446         pr_indent(fp, indent);
3447         fprintf(fp, "lambda = %e\n", sh->lambda);
3448         pr_indent(fp, indent);
3449         fprintf(fp, "buffer size = %" PRId64 "\n", sh->sizeOfTprBody);
3450     }
3451 }