src/gromacs/fileio/tpxio.c

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