src/gromacs/gmxlib/tpxio.c

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