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