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