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