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