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