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