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