src/gmxlib/tpxio.c

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