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