src/gromacs/fileio/tpxio.c

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