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