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