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