src/gromacs/gmxlib/tpxio.c

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