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