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