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