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