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