src/gmxlib/tpxio.c

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