src/gmxlib/tpxio.c

   1 /*
   2  * $Id$
   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 2.0
  11  *
  12  * Copyright (c) 1991-1999
  13  * BIOSON Research Institute, Dept. of Biophysical Chemistry
  14  * University of Groningen, The Netherlands
  15  *
  16  * Please refer to:
  17  * GROMACS: A message-passing parallel molecular dynamics implementation
  18  * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
  19  * Comp. Phys. Comm. 91, 43-56 (1995)
  20  *
  21  * Also check out our WWW page:
  22  * http://md.chem.rug.nl/~gmx
  23  * or e-mail to:
  24  * gromacs@chem.rug.nl
  25  *
  26  * And Hey:
  27  * Green Red Orange Magenta Azure Cyan Skyblue
  28  */
  29 static char *SRCID_tpxio_c = "$Id$";
  30
  31
  32 #include <ctype.h>
  33 #include "sysstuff.h"
  34 #include "smalloc.h"
  35 #include "string2.h"
  36 #include "fatal.h"
  37 #include "macros.h"
  38 #include "assert.h"
  39 #include "names.h"
  40 #include "symtab.h"
  41 #include "futil.h"
  42 #include "filenm.h"
  43 #include "gmxfio.h"
  44 #include "tpxio.h"
  45 #include "confio.h"
  46 #include "atomprop.h"
  47 #include "copyrite.h"
  48
  49 /* This number should be increased whenever the file format changes! */
  50 static int tpx_version = 11;
  51 /* This number should be the most recent incompatible version */
  52 static int tpx_incompatible_version = 9;
  53 /* This is the version of the file we are reading */
  54 static int file_version = 0;
  55
  56
  57 void _do_section(int fp,int key,bool bRead,char *src,int line)
  58 {
  59   char buf[STRLEN];
  60   bool bDbg;
  61
  62   if (fio_getftp(fp) == efTPA) {
  63     if (!bRead) {
  64       do_string(itemstr[key]);
  65       bDbg       = fio_getdebug(fp);
  66       fio_setdebug(fp,FALSE);
  67       do_string(comment_str[key]);
  68       fio_setdebug(fp,bDbg);
  69     }
  70     else {
  71       if (fio_getdebug(fp))
  72         fprintf(stderr,"Looking for section %s (%s, %d)",
  73                 itemstr[key],src,line);
  74
  75       do {
  76         do_string(buf);
  77       } while ((strcasecmp(buf,itemstr[key]) != 0));
  78
  79       if (strcasecmp(buf,itemstr[key]) != 0)
  80         fatal_error(0,"\nCould not find section heading %s",itemstr[key]);
  81       else if (fio_getdebug(fp))
  82         fprintf(stderr," and found it\n");
  83     }
  84   }
  85 }
  86
  87 #define do_section(key,bRead) _do_section(fp,key,bRead,__FILE__,__LINE__)
  88
  89 /**************************************************************
  90  *
  91  * Now the higer level routines that do io of the structures and arrays
  92  *
  93  **************************************************************/
  94 static void do_inputrec(t_inputrec *ir,bool bRead)
  95 {
  96   int  i,j;
  97   bool bDum=TRUE;
  98
  99   if (file_version != tpx_version) {
 100     /* Give a warning about features that are not accessible */
 101     fprintf(stderr,"Note: tpx file_version %d, software version %d\n",
 102             file_version,tpx_version);
 103   }
 104
 105   if (file_version >= 1) {
 106     /* Basic inputrec stuff */
 107     do_int(ir->eI);
 108     do_int(ir->nsteps);
 109     do_int(ir->eBox);
 110     do_int(ir->ns_type);
 111     do_int(ir->nstlist);
 112     do_int(ir->ndelta);
 113     do_int(ir->bDomDecomp);
 114     do_int(ir->decomp_dir);
 115     do_int(ir->nstcomm);
 116     do_int(ir->nstcgsteep);
 117     do_int(ir->nstlog);
 118     do_int(ir->nstxout);
 119     do_int(ir->nstvout);
 120     do_int(ir->nstfout);
 121     do_int(ir->nstenergy);
 122     do_int(ir->nstxtcout);
 123     do_real(ir->init_t);
 124     do_real(ir->delta_t);
 125     do_real(ir->xtcprec);
 126     do_int(ir->solvent_opt);
 127     do_int(ir->nsatoms);
 128     do_int(ir->eBox);
 129     do_real(ir->rlist);
 130     do_int(ir->coulombtype);
 131     do_real(ir->rcoulomb_switch);
 132     do_real(ir->rcoulomb);
 133     do_int(ir->vdwtype);
 134     do_real(ir->rvdw_switch);
 135     do_real(ir->rvdw);
 136     do_int(ir->bDispCorr);
 137     do_real(ir->epsilon_r);
 138     do_int(ir->nkx);
 139     do_int(ir->nky);
 140     do_int(ir->nkz);
 141     do_int(ir->pme_order);
 142     do_real(ir->ewald_rtol);
 143     do_int(ir->bOptFFT);
 144     do_int(ir->bUncStart);
 145     do_int(ir->btc);
 146     do_int(ir->ntcmemory);
 147     do_int(ir->epc);
 148     do_int(ir->npcmemory);
 149     do_real(ir->tau_p);
 150     do_rvec(ir->ref_p);
 151     do_rvec(ir->compress);
 152     do_int(ir->bSimAnn);
 153     do_real(ir->zero_temp_time);
 154     do_real(ir->epsilon_r);
 155     do_real(ir->shake_tol);
 156     do_real(ir->fudgeQQ);
 157     do_int(ir->bPert);
 158     do_real(ir->init_lambda);
 159     do_real(ir->delta_lambda);
 160     do_int(ir->eDisreWeighting);
 161     do_int(ir->bDisreMixed);
 162     do_real(ir->dr_fc);
 163     do_real(ir->dr_tau);
 164     do_int(ir->nstdisreout);
 165     do_real(ir->em_stepsize);
 166     do_real(ir->em_tol);
 167     if (file_version >= 11)
 168       do_int(ir->niter);
 169     else if (bRead) {
 170       ir->niter = 25;
 171       fprintf(stderr,"Note: niter not in run input file, setting it to %d\n",
 172               ir->niter);
 173     }
 174     do_int(ir->eConstrAlg);
 175     do_int(ir->nProjOrder);
 176     do_real(ir->LincsWarnAngle);
 177     do_int(ir->nstLincsout);
 178     do_real(ir->ld_temp);
 179     do_real(ir->ld_fric);
 180     do_int(ir->ld_seed);
 181     do_int(ir->userint1);
 182     do_int(ir->userint2);
 183     do_int(ir->userint3);
 184     do_int(ir->userint4);
 185     do_real(ir->userreal1);
 186     do_real(ir->userreal2);
 187     do_real(ir->userreal3);
 188     do_real(ir->userreal4);
 189
 190     /* grpopts stuff */
 191     do_int(ir->opts.ngtc);
 192     do_int(ir->opts.ngacc);
 193     do_int(ir->opts.ngfrz);
 194     do_int(ir->opts.ngener);
 195     if (bRead) {
 196       snew(ir->opts.nrdf,   ir->opts.ngtc);
 197       snew(ir->opts.ref_t,  ir->opts.ngtc);
 198       snew(ir->opts.tau_t,  ir->opts.ngtc);
 199       snew(ir->opts.nFreeze,ir->opts.ngfrz);
 200       snew(ir->opts.acc,    ir->opts.ngacc);
 201     }
 202     if (ir->opts.ngtc > 0) {
 203       ndo_int (ir->opts.nrdf, ir->opts.ngtc,bDum);
 204       ndo_real(ir->opts.ref_t,ir->opts.ngtc,bDum);
 205       ndo_real(ir->opts.tau_t,ir->opts.ngtc,bDum);
 206     }
 207     if (ir->opts.ngfrz > 0)
 208       ndo_ivec(ir->opts.nFreeze,ir->opts.ngfrz,bDum);
 209     if (ir->opts.ngacc > 0)
 210       ndo_rvec(ir->opts.acc,ir->opts.ngacc);
 211
 212     /* Cosine stuff for electric fields */
 213     for(j=0; (j<DIM); j++) {
 214       do_int  (ir->ex[j].n);
 215       do_int  (ir->et[j].n);
 216       if (bRead) {
 217         snew(ir->ex[j].a,  ir->ex[j].n);
 218         snew(ir->ex[j].phi,ir->ex[j].n);
 219         snew(ir->et[j].a,  ir->et[j].n);
 220         snew(ir->et[j].phi,ir->et[j].n);
 221       }
 222       ndo_real(ir->ex[j].a,  ir->ex[j].n,bDum);
 223       ndo_real(ir->ex[j].phi,ir->ex[j].n,bDum);
 224       ndo_real(ir->et[j].a,  ir->et[j].n,bDum);
 225       ndo_real(ir->et[j].phi,ir->et[j].n,bDum);
 226     }
 227   }
 228   /* set things which are in tpx_version but not in a previous version */
 229
 230   /*
 231   if (file_version < tpx_version) {
 232   } else {
 233   }
 234   */
 235 }
 236
 237 static void do_harm(t_iparams *iparams,bool bRead)
 238 {
 239   do_real(iparams->harmonic.rA);
 240   do_real(iparams->harmonic.krA);
 241   do_real(iparams->harmonic.rB);
 242   do_real(iparams->harmonic.krB);
 243 }
 244
 245 void do_iparams(t_functype ftype,t_iparams *iparams,bool bRead)
 246 {
 247   int i;
 248   bool bDum;
 249
 250   if (!bRead)
 251     set_comment(interaction_function[ftype].name);
 252   switch (ftype) {
 253   case F_ANGLES:
 254   case F_G96ANGLES:
 255   case F_BONDS:
 256   case F_G96BONDS:
 257   case F_IDIHS:
 258   case F_ANGRES:
 259   case F_ANGRESZ:
 260     do_harm(iparams,bRead);
 261     break;
 262   case F_BHAM:
 263     do_real(iparams->bham.a);
 264     do_real(iparams->bham.b);
 265     do_real(iparams->bham.c);
 266     break;
 267   case F_MORSE:
 268     do_real(iparams->morse.b0);
 269     do_real(iparams->morse.cb);
 270     do_real(iparams->morse.beta);
 271     break;
 272   case F_WPOL:
 273     do_real(iparams->wpol.kx);
 274     do_real(iparams->wpol.ky);
 275     do_real(iparams->wpol.kz);
 276     do_real(iparams->wpol.rOH);
 277     do_real(iparams->wpol.rHH);
 278     do_real(iparams->wpol.rOD);
 279     break;
 280   case F_LJ:
 281     do_real(iparams->lj.c6);
 282     do_real(iparams->lj.c12);
 283     break;
 284   case F_LJ14:
 285     do_real(iparams->lj14.c6A);
 286     do_real(iparams->lj14.c12A);
 287     do_real(iparams->lj14.c6B);
 288     do_real(iparams->lj14.c12B);
 289     break;
 290   case F_PDIHS:
 291     do_real(iparams->pdihs.phiA);
 292     do_real(iparams->pdihs.cpA);
 293     do_real(iparams->pdihs.phiB);
 294     do_real(iparams->pdihs.cpB);
 295     do_int (iparams->pdihs.mult);
 296     break;
 297   case F_DISRES:
 298     do_int (iparams->disres.index);
 299     do_int (iparams->disres.type);
 300     do_real(iparams->disres.low);
 301     do_real(iparams->disres.up1);
 302     do_real(iparams->disres.up2);
 303     do_real(iparams->disres.fac);
 304     break;
 305   case F_POSRES:
 306     do_rvec(iparams->posres.pos0);
 307     do_rvec(iparams->posres.fc);
 308     break;
 309   case F_RBDIHS:
 310     ndo_real(iparams->rbdihs.rbc,NR_RBDIHS,bDum);
 311     break;
 312   case F_SHAKE:
 313   case F_SHAKENC:
 314     do_real(iparams->shake.dA);
 315     do_real(iparams->shake.dB);
 316     break;
 317   case F_SETTLE:
 318     do_real(iparams->settle.doh);
 319     do_real(iparams->settle.dhh);
 320     break;
 321   case F_DUMMY2:
 322     do_real(iparams->dummy.a);
 323     break;
 324   case F_DUMMY3:
 325   case F_DUMMY3FD:
 326   case F_DUMMY3FAD:
 327     do_real(iparams->dummy.a);
 328     do_real(iparams->dummy.b);
 329     break;
 330   case F_DUMMY3OUT:
 331   case F_DUMMY4FD:
 332     do_real(iparams->dummy.a);
 333     do_real(iparams->dummy.b);
 334     do_real(iparams->dummy.c);
 335     break;
 336   default:
 337     fatal_error(0,"unknown function type %d (%s) in %s line %d",
 338                 ftype,interaction_function[ftype].name,__FILE__,__LINE__);
 339   }
 340   if (!bRead)
 341     unset_comment();
 342 }
 343
 344 static void do_ilist(t_ilist *ilist,bool bRead,char *name)
 345 {
 346   int i;
 347   bool bDum=TRUE;
 348
 349   if (!bRead)
 350     set_comment(name);
 351   ndo_int(ilist->multinr,MAXPROC,bDum);
 352   do_int (ilist->nr);
 353   if (bRead)
 354     snew(ilist->iatoms,ilist->nr);
 355   ndo_int(ilist->iatoms,ilist->nr,bDum);
 356   if (!bRead)
 357     unset_comment();
 358 }
 359
 360 static void do_idef(t_idef *idef,bool bRead)
 361 {
 362   int i,j;
 363   bool bDum=TRUE;
 364
 365   do_int(idef->atnr);
 366   do_int(idef->pid);
 367   do_int(idef->ntypes);
 368   if (bRead) {
 369     snew(idef->functype,idef->ntypes);
 370     snew(idef->iparams,idef->ntypes);
 371   }
 372   ndo_int(idef->functype,idef->ntypes,bDum);
 373
 374   for (i=0; (i<idef->ntypes); i++)
 375     do_iparams(idef->functype[i],&idef->iparams[i],bRead);
 376
 377   for(j=0; (j<F_NRE); j++) {
 378     if ((bRead && (file_version < 6)) &&
 379         ((j == F_G96ANGLES) || (j == F_G96BONDS))) {
 380       fprintf(stderr,"Warning: file_version %d < 6: no GROMOS96 Force field\n"
 381               ,file_version);
 382       idef->il[j].nr         = 0;
 383       idef->il[j].multinr[0] = 0;
 384       idef->il[j].iatoms     = NULL;
 385     }
 386     else
 387       do_ilist(&idef->il[j],bRead,interaction_function[j].name);
 388   }
 389 }
 390
 391 static void do_block(t_block *block,bool bRead)
 392 {
 393   int i;
 394   bool bDum=TRUE;
 395
 396   ndo_int(block->multinr,MAXPROC,bDum);
 397   do_int (block->nr);
 398   do_int (block->nra);
 399   if (bRead) {
 400     snew(block->index,block->nr+1);
 401     snew(block->a,block->nra);
 402   }
 403   ndo_int(block->index,block->nr+1,bDum);
 404   ndo_int(block->a,block->nra,bDum);
 405 }
 406
 407 static void do_atom(t_atom *atom,bool bRead)
 408 {
 409   int ngrp=egcNR;
 410
 411   do_real (atom->m);
 412   do_real (atom->q);
 413   do_real (atom->mB);
 414   do_real (atom->qB);
 415   do_ushort(atom->type);
 416   do_ushort(atom->typeB);
 417   do_int (atom->ptype);
 418   do_int (atom->resnr);
 419   do_nuchar(atom->grpnr,ngrp);
 420 }
 421
 422 static void do_grps(int ngrp,t_grps grps[],bool bRead)
 423 {
 424   int i,j;
 425   bool bDum=TRUE;
 426
 427   for(j=0; (j<ngrp); j++) {
 428     do_int (grps[j].nr);
 429     if (bRead)
 430       snew(grps[j].nm_ind,grps[j].nr);
 431     ndo_int(grps[j].nm_ind,grps[j].nr,bDum);
 432   }
 433 }
 434
 435 static void do_symstr(char ***nm,bool bRead,t_symtab *symtab)
 436 {
 437   int ls;
 438
 439   if (bRead) {
 440     do_int(ls);
 441     *nm = get_symtab_handle(symtab,ls);
 442   }
 443   else {
 444     ls = lookup_symtab(symtab,*nm);
 445     do_int(ls);
 446   }
 447 }
 448
 449 static void do_strstr(int nstr,char ***nm,bool bRead,t_symtab *symtab)
 450 {
 451   int  j;
 452
 453   for (j=0; (j<nstr); j++)
 454     do_symstr(&(nm[j]),bRead,symtab);
 455 }
 456
 457 static void do_atoms(t_atoms *atoms,bool bRead,t_symtab *symtab)
 458 {
 459   int i;
 460
 461   do_int(atoms->nr);
 462   do_int(atoms->nres);
 463   do_int(atoms->ngrpname);
 464   if (bRead) {
 465     snew(atoms->atom,atoms->nr);
 466     snew(atoms->atomname,atoms->nr);
 467     snew(atoms->resname,atoms->nres);
 468     snew(atoms->grpname,atoms->ngrpname);
 469     atoms->pdbinfo = NULL;
 470   }
 471   for(i=0; (i<atoms->nr); i++)
 472     do_atom(&atoms->atom[i],bRead);
 473   do_strstr(atoms->nr,atoms->atomname,bRead,symtab);
 474   do_strstr(atoms->nres,atoms->resname,bRead,symtab);
 475   do_strstr(atoms->ngrpname,atoms->grpname,bRead,symtab);
 476
 477   do_grps(egcNR,atoms->grps,bRead);
 478
 479   do_block(&(atoms->excl),bRead);
 480 }
 481
 482 static void do_symtab(t_symtab *symtab,bool bRead)
 483 {
 484   int i,nr;
 485   t_symbuf *symbuf;
 486   char buf[STRLEN];
 487
 488   do_int(symtab->nr);
 489   nr     = symtab->nr;
 490   if (bRead) {
 491     snew(symtab->symbuf,1);
 492     symbuf = symtab->symbuf;
 493     symbuf->bufsize = nr;
 494     snew(symbuf->buf,nr);
 495     for (i=0; (i<nr); i++) {
 496       do_string(buf);
 497       symbuf->buf[i]=strdup(buf);
 498     }
 499   }
 500   else {
 501     symbuf = symtab->symbuf;
 502     while (symbuf!=NULL) {
 503       for (i=0; (i<symbuf->bufsize) && (i<nr); i++)
 504         do_string(symbuf->buf[i]);
 505       nr-=i;
 506       symbuf=symbuf->next;
 507     }
 508     if (nr != 0)
 509       fatal_error(0,"nr of symtab strings left: %d",nr);
 510   }
 511 }
 512
 513 static void make_chain_identifiers(t_atoms *atoms,t_block *mols)
 514 {
 515   int m,a,a0,a1;
 516   char c,chain;
 517 #define CHAIN_MIN_ATOMS 15
 518
 519   chain='A';
 520   for(m=0; m<mols->nr; m++) {
 521     a0=mols->index[m];
 522     a1=mols->index[m+1];
 523     if ((a1-a0 >= CHAIN_MIN_ATOMS) && (chain <= 'Z')) {
 524       c=chain;
 525       chain++;
 526     } else
 527       c=' ';
 528     for(a=a0; a<a1; a++)
 529       atoms->atom[mols->a[a]].chain=c;
 530   }
 531   if (chain == 'B')
 532     for(a=0; a<atoms->nr; a++)
 533       atoms->atom[a].chain=' ';
 534 }
 535
 536 static void do_top(t_topology *top,bool bRead)
 537 {
 538   int  i;
 539
 540   do_symtab(&(top->symtab),bRead);
 541   do_symstr(&(top->name),bRead,&(top->symtab));
 542   do_atoms (&(top->atoms),bRead,&(top->symtab));
 543   do_idef  (&(top->idef),bRead);
 544   for(i=0; (i<ebNR); i++)
 545     do_block(&(top->blocks[i]),bRead);
 546   if (bRead) {
 547     close_symtab(&(top->symtab));
 548     make_chain_identifiers(&(top->atoms),&(top->blocks[ebMOLS]));
 549   }
 550 }
 551
 552 static void do_tpxheader(int fp,bool bRead,t_tpxheader *tpx)
 553 {
 554   char  buf[STRLEN];
 555   bool  bDouble;
 556   int   precision;
 557
 558   fio_select(fp);
 559   fio_setdebug(fp,bDebugMode());
 560
 561   /* NEW! XDR tpb file */
 562   precision = sizeof(real);
 563   if (bRead) {
 564     do_string(buf);
 565     if (strncmp(buf,"VERSION",7))
 566       fatal_error(0,"Can not read file %s,\n"
 567                   "             this file is from a Gromacs version which is older than 2.0\n"
 568                   "             Make a new one with grompp or use a gro or pdb file, if possible",
 569                   fio_getname(fp));
 570     do_int(precision);
 571     bDouble = (precision == sizeof(double));
 572     if ((precision != sizeof(float)) && !bDouble)
 573       fatal_error(0,"Unknown precision in file %s: real is %d bytes "
 574                   "instead of %d or %d",
 575                   fio_getname(fp),precision,sizeof(float),sizeof(double));
 576     fio_setprecision(fp,bDouble);
 577     fprintf(stderr,"Reading file %s, %s (%s precision)\n",
 578             fio_getname(fp),buf,bDouble ? "double" : "single");
 579   }
 580   else {
 581     do_string(GromacsVersion());
 582     bDouble = (precision == sizeof(double));
 583     fio_setprecision(fp,bDouble);
 584     do_int(precision);
 585     file_version = tpx_version;
 586   }
 587
 588   /* Check versions! */
 589   do_int(file_version);
 590   if ((file_version <= tpx_incompatible_version) ||
 591       (file_version > tpx_version))
 592     fatal_error(0,"reading tpx file (%s) version %d with version %d program",
 593                 fio_getname(fp),file_version,tpx_version);
 594
 595   do_section(eitemHEADER,bRead);
 596   do_int (tpx->natoms);
 597   do_int (tpx->step);
 598   do_real(tpx->t);
 599   do_real(tpx->lambda);
 600   do_int (tpx->bIr);
 601   do_int (tpx->bTop);
 602   do_int (tpx->bX);
 603   do_int (tpx->bV);
 604   do_int (tpx->bF);
 605   do_int (tpx->bBox);
 606 }
 607
 608 static void do_tpx(int fp,bool bRead,int *step,real *t,real *lambda,
 609                    t_inputrec *ir,rvec *box,int *natoms,
 610                    rvec *x,rvec *v,rvec *f,t_topology *top)
 611 {
 612   t_tpxheader tpx;
 613   t_inputrec  dum_ir;
 614   t_topology  dum_top;
 615
 616   if (!bRead) {
 617     tpx.natoms = *natoms;
 618     tpx.step   = *step;
 619     tpx.t      = *t;
 620     tpx.lambda = *lambda;
 621     tpx.bIr  = (ir  != NULL);
 622     tpx.bTop = (top != NULL);
 623     tpx.bX   = (x   != NULL);
 624     tpx.bV   = (v   != NULL);
 625     tpx.bF   = (f   != NULL);
 626     tpx.bBox = (box != NULL);
 627   }
 628
 629   do_tpxheader(fp,bRead,&tpx);
 630
 631   if (bRead) {
 632     *natoms = tpx.natoms;
 633     *step   = tpx.step;
 634     *t      = tpx.t;
 635     *lambda = tpx.lambda;
 636   }
 637
 638 #define do_test(b,p) if (bRead && (p!=NULL) && !b) fatal_error(0,"No %s in %s",#p,fio_getname(fp))
 639
 640   do_test(tpx.bBox,box);
 641   do_section(eitemBOX,bRead);
 642   if (tpx.bBox) ndo_rvec(box,DIM);
 643
 644   do_test(tpx.bIr,ir);
 645   do_section(eitemIR,bRead);
 646   if (tpx.bIr) {
 647     if (ir)
 648       do_inputrec(ir,bRead);
 649     else {
 650       init_inputrec(&dum_ir);
 651       do_inputrec  (&dum_ir,bRead);
 652       done_inputrec(&dum_ir);
 653     }
 654   }
 655   do_test(tpx.bTop,top);
 656   do_section(eitemTOP,bRead);
 657   if (tpx.bTop) {
 658     if (top)
 659       do_top(top,bRead);
 660     else {
 661       do_top(&dum_top,bRead);
 662       done_top(&dum_top);
 663     }
 664   }
 665   do_test(tpx.bX,x);
 666   do_section(eitemX,bRead);
 667   if (tpx.bX) ndo_rvec(x,*natoms);
 668
 669   do_test(tpx.bV,v);
 670   do_section(eitemV,bRead);
 671   if (tpx.bV) ndo_rvec(v,*natoms);
 672
 673   do_test(tpx.bF,f);
 674   do_section(eitemF,bRead);
 675   if (tpx.bF) ndo_rvec(f,*natoms);
 676 }
 677
 678 /************************************************************
 679  *
 680  *  The following routines are the exported ones
 681  *
 682  ************************************************************/
 683
 684 int open_tpx(char *fn,char *mode)
 685 {
 686   return fio_open(fn,mode);
 687 }
 688
 689 void close_tpx(int fp)
 690 {
 691   fio_close(fp);
 692 }
 693
 694 void read_tpxheader(char *fn,t_tpxheader *tpx)
 695 {
 696   int fp;
 697
 698   fp = open_tpx(fn,"r");
 699   do_tpxheader(fp,TRUE,tpx);
 700   close_tpx(fp);
 701 }
 702
 703 void write_tpx(char *fn,int step,real t,real lambda,
 704                t_inputrec *ir,rvec *box,int natoms,
 705                rvec *x,rvec *v,rvec *f,t_topology *top)
 706 {
 707   int fp;
 708
 709   fp = open_tpx(fn,"w");
 710   do_tpx(fp,FALSE,&step,&t,&lambda,ir,box,&natoms,x,v,f,top);
 711   close_tpx(fp);
 712 }
 713
 714 void read_tpx(char *fn,int *step,real *t,real *lambda,
 715               t_inputrec *ir,rvec *box,int *natoms,
 716               rvec *x,rvec *v,rvec *f,t_topology *top)
 717 {
 718   int fp;
 719
 720   fp = open_tpx(fn,"r");
 721   do_tpx(fp,TRUE,step,t,lambda,ir,box,natoms,x,v,f,top);
 722   close_tpx(fp);
 723 }
 724
 725 void fwrite_tpx(int fp,int step,real t,real lambda,
 726                 t_inputrec *ir,rvec *box,int natoms,
 727                 rvec *x,rvec *v,rvec *f,t_topology *top)
 728 {
 729   do_tpx(fp,FALSE,&step,&t,&lambda,ir,box,&natoms,x,v,f,top);
 730 }
 731
 732
 733 void fread_tpx(int fp,int *step,real *t,real *lambda,
 734                t_inputrec *ir,rvec *box,int *natoms,
 735                rvec *x,rvec *v,rvec *f,t_topology *top)
 736 {
 737   do_tpx(fp,TRUE,step,t,lambda,ir,box,natoms,x,v,f,top);
 738 }
 739
 740 bool fn2bTPX(char *file)
 741 {
 742   switch (fn2ftp(file)) {
 743   case efTPR:
 744   case efTPB:
 745   case efTPA:
 746     return TRUE;
 747   default:
 748     return FALSE;
 749   }
 750 }
 751
 752 bool read_tps_conf(char *infile,char *title,t_topology *top,
 753                    rvec **x,rvec **v,matrix box,bool bMass)
 754 {
 755   t_tpxheader  header;
 756   real         t,lambda;
 757   int          natoms,step,i;
 758   bool         bTop;
 759
 760   bTop=fn2bTPX(infile);
 761   if (bTop) {
 762     read_tpxheader(infile,&header);
 763     snew(*x,header.natoms);
 764     if (v)
 765       snew(*v,header.natoms);
 766     read_tpx(infile,&step,&t,&lambda,NULL,box,&natoms,
 767              *x,(v==NULL) ? NULL : *v,NULL,top);
 768     strcpy(title,*top->name);
 769   }
 770   else {
 771     get_stx_coordnum(infile,&natoms);
 772     init_t_atoms(&top->atoms,natoms,FALSE);
 773     snew(*x,natoms);
 774     if (v)
 775       snew(*v,natoms);
 776     read_stx_conf(infile,title,&top->atoms,*x,(v==NULL) ? NULL : *v,box);
 777     if (bMass)
 778       for(i=0; i<natoms; i++)
 779         top->atoms.atom[i].m =
 780           get_mass(*top->atoms.resname[top->atoms.atom[i].resnr],
 781                    *top->atoms.atomname[i]);
 782     top->idef.ntypes=-1;
 783   }
 784
 785   return bTop;
 786 }