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