src/gmxlib/tpxio.c

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