src/gromacs/gmxlib/confio.c

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  *                        VERSION 3.2.0
  10  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  11  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  12  * Copyright (c) 2001-2004, The GROMACS development team,
  13  * check out http://www.gromacs.org for more information.
  14
  15  * This program is free software; you can redistribute it and/or
  16  * modify it under the terms of the GNU General Public License
  17  * as published by the Free Software Foundation; either version 2
  18  * of the License, or (at your option) any later version.
  19  *
  20  * If you want to redistribute modifications, please consider that
  21  * scientific software is very special. Version control is crucial -
  22  * bugs must be traceable. We will be happy to consider code for
  23  * inclusion in the official distribution, but derived work must not
  24  * be called official GROMACS. Details are found in the README & COPYING
  25  * files - if they are missing, get the official version at www.gromacs.org.
  26  *
  27  * To help us fund GROMACS development, we humbly ask that you cite
  28  * the papers on the package - you can find them in the top README file.
  29  *
  30  * For more info, check our website at http://www.gromacs.org
  31  *
  32  * And Hey:
  33  * GROningen Mixture of Alchemy and Childrens' Stories
  34  */
  35 #ifdef HAVE_CONFIG_H
  36 #include <config.h>
  37 #endif
  38
  39 #include <math.h>
  40 #include "sysstuff.h"
  41 #include "typedefs.h"
  42 #include "smalloc.h"
  43 #include "sysstuff.h"
  44 #include <errno.h>
  45 #include "macros.h"
  46 #include "string2.h"
  47 #include "confio.h"
  48 #include "vec.h"
  49 #include "symtab.h"
  50 #include "futil.h"
  51 #include "xdrf.h"
  52 #include "filenm.h"
  53 #include "pdbio.h"
  54 #include "tpxio.h"
  55 #include "gmx_fatal.h"
  56 #include "copyrite.h"
  57 #include "filenm.h"
  58 #include "statutil.h"
  59 #include "pbc.h"
  60 #include "mtop_util.h"
  61 #include "gmxfio.h"
  62
  63 #define CHAR_SHIFT 24
  64
  65 static int read_g96_pos(char line[],t_symtab *symtab,
  66                         FILE *fp,const char *infile,
  67                         t_trxframe *fr)
  68 {
  69   t_atoms *atoms;
  70   gmx_bool   bEnd;
  71   int    nwanted,natoms,atnr,resnr=0,oldres,newres,shift;
  72   char   anm[STRLEN],resnm[STRLEN];
  73   char   c1,c2;
  74   double db1,db2,db3;
  75
  76   nwanted = fr->natoms;
  77
  78   atoms = fr->atoms;
  79
  80   natoms = 0;
  81
  82   if (fr->bX) {
  83     if (fr->bAtoms)
  84       shift = CHAR_SHIFT;
  85     else
  86       shift = 0;
  87     newres  = -1;
  88     oldres  = -666; /* Unlikely number for the first residue! */
  89     bEnd    = FALSE;
  90     while (!bEnd && fgets2(line,STRLEN,fp)) {
  91       bEnd = (strncmp(line,"END",3) == 0);
  92       if (!bEnd  && (line[0] != '#')) {
  93         if (sscanf(line+shift,"%15lf%15lf%15lf",&db1,&db2,&db3) != 3)
  94           gmx_fatal(FARGS,"Did not find 3 coordinates for atom %d in %s\n",
  95                       natoms+1,infile);
  96         if ((nwanted != -1) && (natoms >= nwanted))
  97           gmx_fatal(FARGS,
  98                       "Found more coordinates (%d) in %s than expected %d\n",
  99                       natoms,infile,nwanted);
 100         if (atoms) {
 101           if (fr->bAtoms &&
 102               (sscanf(line,"%5d%c%5s%c%5s%7d",&resnr,&c1,resnm,&c2,anm,&atnr)
 103                != 6)) {
 104             if (oldres>=0)
 105               resnr = oldres;
 106             else {
 107               resnr    = 1;
 108               strncpy(resnm,"???",sizeof(resnm)-1);
 109             }
 110             strncpy(anm,"???",sizeof(anm)-1);
 111           }
 112           atoms->atomname[natoms]=put_symtab(symtab,anm);
 113           if (resnr != oldres) {
 114             oldres = resnr;
 115             newres++;
 116             if (newres >= atoms->nr)
 117               gmx_fatal(FARGS,"More residues than atoms in %s (natoms = %d)",
 118                           infile,atoms->nr);
 119             atoms->atom[natoms].resind = newres;
 120             if (newres+1 > atoms->nres) {
 121               atoms->nres = newres+1;
 122             }
 123             t_atoms_set_resinfo(atoms,natoms,symtab,resnm,resnr,' ',0,' ');
 124           } else {
 125             atoms->atom[natoms].resind = newres;
 126           }
 127         }
 128         if (fr->x) {
 129           fr->x[natoms][0] = db1;
 130           fr->x[natoms][1] = db2;
 131           fr->x[natoms][2] = db3;
 132         }
 133         natoms++;
 134       }
 135     }
 136     if ((nwanted != -1) && natoms != nwanted)
 137       fprintf(stderr,
 138               "Warning: found less coordinates (%d) in %s than expected %d\n",
 139               natoms,infile,nwanted);
 140   }
 141
 142   fr->natoms = natoms;
 143
 144   return natoms;
 145 }
 146
 147 static int read_g96_vel(char line[],FILE *fp,const char *infile,
 148                         t_trxframe *fr)
 149 {
 150   gmx_bool   bEnd;
 151   int    nwanted,natoms=-1,shift;
 152   double db1,db2,db3;
 153
 154   nwanted = fr->natoms;
 155
 156   if (fr->v && fr->bV) {
 157     if (strcmp(line,"VELOCITYRED") == 0)
 158       shift = 0;
 159     else
 160       shift = CHAR_SHIFT;
 161     natoms = 0;
 162     bEnd = FALSE;
 163     while (!bEnd && fgets2(line,STRLEN,fp)) {
 164       bEnd = (strncmp(line,"END",3) == 0);
 165       if (!bEnd && (line[0] != '#')) {
 166         if (sscanf(line+shift,"%15lf%15lf%15lf",&db1,&db2,&db3) != 3)
 167           gmx_fatal(FARGS,"Did not find 3 velocities for atom %d in %s",
 168                       natoms+1,infile);
 169         if ((nwanted != -1) && (natoms >= nwanted))
 170           gmx_fatal(FARGS,"Found more velocities (%d) in %s than expected %d\n",
 171                       natoms,infile,nwanted);
 172         if (fr->v) {
 173           fr->v[natoms][0] = db1;
 174           fr->v[natoms][1] = db2;
 175           fr->v[natoms][2] = db3;
 176         }
 177         natoms++;
 178       }
 179     }
 180     if ((nwanted != -1) && (natoms != nwanted))
 181       fprintf(stderr,
 182               "Warning: found less velocities (%d) in %s than expected %d\n",
 183               natoms,infile,nwanted);
 184   }
 185
 186   return natoms;
 187 }
 188
 189 int read_g96_conf(FILE *fp,const char *infile,t_trxframe *fr, char *line)
 190 {
 191   t_symtab *symtab=NULL;
 192   gmx_bool   bAtStart,bTime,bAtoms,bPos,bVel,bBox,bEnd,bFinished;
 193   int    natoms,nbp;
 194   double db1,db2,db3,db4,db5,db6,db7,db8,db9;
 195
 196   bAtStart = (ftell(fp) == 0);
 197
 198   clear_trxframe(fr,FALSE);
 199
 200   if (!symtab) {
 201     snew(symtab,1);
 202     open_symtab(symtab);
 203   }
 204
 205   natoms=0;
 206
 207   if (bAtStart) {
 208     while ( !fr->bTitle && fgets2(line,STRLEN,fp))
 209       fr->bTitle = (strcmp(line,"TITLE") == 0);
 210     if (fr->title == NULL) {
 211       fgets2(line,STRLEN,fp);
 212       fr->title = strdup(line);
 213     }
 214     bEnd = FALSE;
 215     while (!bEnd && fgets2(line,STRLEN,fp))
 216       bEnd = (strcmp(line,"END") == 0);
 217     fgets2(line,STRLEN,fp);
 218   }
 219
 220   /* Do not get a line if we are not at the start of the file, *
 221    * because without a parameter file we don't know what is in *
 222    * the trajectory and we have already read the line in the   *
 223    * previous call (VERY DIRTY).                               */
 224   bFinished = FALSE;
 225   do {
 226     bTime  = (strcmp(line,"TIMESTEP") == 0);
 227     bAtoms = (strcmp(line,"POSITION") == 0);
 228     bPos   = (bAtoms || (strcmp(line,"POSITIONRED") == 0));
 229     bVel   = (strncmp(line,"VELOCITY",8) == 0);
 230     bBox   = (strcmp(line,"BOX") == 0);
 231     if (bTime) {
 232       if (!fr->bTime && !fr->bX) {
 233         fr->bStep = bTime;
 234         fr->bTime = bTime;
 235         do
 236           bFinished = (fgets2(line,STRLEN,fp) == NULL);
 237         while (!bFinished && (line[0] == '#'));
 238         sscanf(line,"%15d%15lf",&(fr->step),&db1);
 239         fr->time = db1;
 240       } else
 241         bFinished = TRUE;
 242     }
 243     if (bPos) {
 244       if (!fr->bX) {
 245         fr->bAtoms = bAtoms;
 246         fr->bX     = bPos;
 247         natoms = read_g96_pos(line,symtab,fp,infile,fr);
 248       } else
 249         bFinished = TRUE;
 250     }
 251     if (fr->v && bVel) {
 252       fr->bV = bVel;
 253       natoms = read_g96_vel(line,fp,infile,fr);
 254     }
 255     if (bBox) {
 256       fr->bBox = bBox;
 257       clear_mat(fr->box);
 258       bEnd = FALSE;
 259       while (!bEnd && fgets2(line,STRLEN,fp)) {
 260         bEnd = (strncmp(line,"END",3) == 0);
 261         if (!bEnd && (line[0] != '#')) {
 262           nbp = sscanf(line,"%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf",
 263                        &db1,&db2,&db3,&db4,&db5,&db6,&db7,&db8,&db9);
 264           if (nbp < 3)
 265             gmx_fatal(FARGS,"Found a BOX line, but no box in %s",infile);
 266           fr->box[XX][XX] = db1;
 267           fr->box[YY][YY] = db2;
 268           fr->box[ZZ][ZZ] = db3;
 269           if (nbp == 9) {
 270             fr->box[XX][YY] = db4;
 271             fr->box[XX][ZZ] = db5;
 272             fr->box[YY][XX] = db6;
 273             fr->box[YY][ZZ] = db7;
 274             fr->box[ZZ][XX] = db8;
 275             fr->box[ZZ][YY] = db9;
 276           }
 277         }
 278       }
 279       bFinished = TRUE;
 280     }
 281   } while (!bFinished && fgets2(line,STRLEN,fp));
 282
 283   free_symtab(symtab);
 284
 285   fr->natoms = natoms;
 286
 287   return natoms;
 288 }
 289
 290 void write_g96_conf(FILE *out,t_trxframe *fr,
 291                     int nindex,const atom_id *index)
 292 {
 293   t_atoms *atoms;
 294   int nout,i,a;
 295
 296   atoms = fr->atoms;
 297
 298   if (index)
 299     nout = nindex;
 300   else
 301     nout = fr->natoms;
 302
 303   if (fr->bTitle)
 304     fprintf(out,"TITLE\n%s\nEND\n",fr->title);
 305   if (fr->bStep || fr->bTime)
 306     /* Officially the time format is %15.9, which is not enough for 10 ns */
 307     fprintf(out,"TIMESTEP\n%15d%15.6f\nEND\n",fr->step,fr->time);
 308   if (fr->bX) {
 309     if (fr->bAtoms) {
 310       fprintf(out,"POSITION\n");
 311       for(i=0; i<nout; i++) {
 312         if (index) a = index[i]; else a = i;
 313         fprintf(out,"%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
 314                 (atoms->resinfo[atoms->atom[a].resind].nr) % 100000,
 315                 *atoms->resinfo[atoms->atom[a].resind].name,
 316                 *atoms->atomname[a],(i+1) % 10000000,
 317                 fr->x[a][XX],fr->x[a][YY],fr->x[a][ZZ]);
 318       }
 319     } else {
 320       fprintf(out,"POSITIONRED\n");
 321       for(i=0; i<nout; i++) {
 322         if (index) a = index[i]; else a = i;
 323         fprintf(out,"%15.9f%15.9f%15.9f\n",
 324                 fr->x[a][XX],fr->x[a][YY],fr->x[a][ZZ]);
 325       }
 326     }
 327     fprintf(out,"END\n");
 328   }
 329   if (fr->bV) {
 330     if (fr->bAtoms) {
 331       fprintf(out,"VELOCITY\n");
 332       for(i=0; i<nout; i++) {
 333         if (index) a = index[i]; else a = i;
 334         fprintf(out,"%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
 335                 (atoms->resinfo[atoms->atom[a].resind].nr) % 100000,
 336                 *atoms->resinfo[atoms->atom[a].resind].name,
 337                 *atoms->atomname[a],(i+1) % 10000000,
 338                 fr->v[a][XX],fr->v[a][YY],fr->v[a][ZZ]);
 339       }
 340     } else {
 341       fprintf(out,"VELOCITYRED\n");
 342       for(i=0; i<nout; i++) {
 343         if (index) a = index[i]; else a = i;
 344         fprintf(out,"%15.9f%15.9f%15.9f\n",
 345                 fr->v[a][XX],fr->v[a][YY],fr->v[a][ZZ]);
 346       }
 347     }
 348     fprintf(out,"END\n");
 349   }
 350   if (fr->bBox) {
 351     fprintf(out,"BOX\n");
 352     fprintf(out,"%15.9f%15.9f%15.9f",
 353             fr->box[XX][XX],fr->box[YY][YY],fr->box[ZZ][ZZ]);
 354     if (fr->box[XX][YY] || fr->box[XX][ZZ] || fr->box[YY][XX] ||
 355         fr->box[YY][ZZ] || fr->box[ZZ][XX] ||fr->box[ZZ][YY])
 356       fprintf(out,"%15.9f%15.9f%15.9f%15.9f%15.9f%15.9f",
 357               fr->box[XX][YY],fr->box[XX][ZZ],fr->box[YY][XX],
 358               fr->box[YY][ZZ],fr->box[ZZ][XX],fr->box[ZZ][YY]);
 359     fprintf(out,"\n");
 360     fprintf(out,"END\n");
 361   }
 362 }
 363
 364 static int get_espresso_word(FILE *fp,char word[])
 365 {
 366   int  ret,nc,i;
 367
 368   ret = 0;
 369   nc = 0;
 370
 371   do {
 372     i = fgetc(fp);
 373     if (i != EOF) {
 374       if (i == ' ' || i == '\n' || i == '\t') {
 375         if (nc > 0)
 376           ret = 1;
 377       } else if (i == '{') {
 378         if (nc == 0)
 379           word[nc++] = '{';
 380         ret = 2;
 381       } else if (i == '}') {
 382         if (nc == 0)
 383           word[nc++] = '}';
 384         ret = 3;
 385       } else {
 386         word[nc++] = (char)i;
 387       }
 388     }
 389   } while (i != EOF && ret == 0);
 390
 391   word[nc] = '\0';
 392
 393   /*  printf("'%s'\n",word); */
 394
 395   return ret;
 396 }
 397
 398 static int check_open_parenthesis(FILE *fp,int r,
 399                                   const char *infile,const char *keyword)
 400 {
 401   int level_inc;
 402   char word[STRLEN];
 403
 404   level_inc = 0;
 405   if (r == 2) {
 406     level_inc++;
 407   } else {
 408     r = get_espresso_word(fp,word);
 409     if (r == 2)
 410       level_inc++;
 411     else
 412       gmx_fatal(FARGS,"Expected '{' after '%s' in file '%s'",
 413                 keyword,infile);
 414   }
 415
 416   return level_inc;
 417 }
 418
 419 static int check_close_parenthesis(FILE *fp,int r,
 420                                   const char *infile,const char *keyword)
 421 {
 422   int level_inc;
 423   char word[STRLEN];
 424
 425   level_inc = 0;
 426   if (r == 3) {
 427     level_inc--;
 428   } else {
 429     r = get_espresso_word(fp,word);
 430     if (r == 3)
 431       level_inc--;
 432     else
 433       gmx_fatal(FARGS,"Expected '}' after section '%s' in file '%s'",
 434                 keyword,infile);
 435   }
 436
 437   return level_inc;
 438 }
 439
 440 enum { espID, espPOS, espTYPE, espQ, espV, espF, espMOLECULE, espNR };
 441 const char *esp_prop[espNR] = { "id", "pos", "type", "q", "v", "f",
 442                                 "molecule" };
 443
 444 static void read_espresso_conf(const char *infile,
 445                                t_atoms *atoms,rvec x[],rvec *v,matrix box)
 446 {
 447   t_symtab *symtab=NULL;
 448   FILE *fp;
 449   char word[STRLEN],buf[STRLEN];
 450   int  natoms,level,npar,r,nprop,p,i,m,molnr;
 451   int  prop[32];
 452   double d;
 453   gmx_bool bFoundParticles,bFoundProp,bFoundVariable,bMol;
 454
 455   if (!symtab) {
 456     snew(symtab,1);
 457     open_symtab(symtab);
 458   }
 459
 460   clear_mat(box);
 461
 462   fp = gmx_fio_fopen(infile,"r");
 463
 464   bFoundParticles = FALSE;
 465   bFoundVariable = FALSE;
 466   bMol = FALSE;
 467   level = 0;
 468   while ((r=get_espresso_word(fp,word))) {
 469     if (level==1 && strcmp(word,"particles")==0 && !bFoundParticles) {
 470       bFoundParticles = TRUE;
 471       level += check_open_parenthesis(fp,r,infile,"particles");
 472       nprop = 0;
 473       while (level == 2 && (r=get_espresso_word(fp,word))) {
 474         bFoundProp = FALSE;
 475         for(p=0; p<espNR; p++) {
 476           if (strcmp(word,esp_prop[p]) == 0) {
 477             bFoundProp = TRUE;
 478             prop[nprop++] = p;
 479             /* printf("  prop[%d] = %s\n",nprop-1,esp_prop[prop[nprop-1]]); */
 480           }
 481         }
 482         if (!bFoundProp && word[0] != '}') {
 483           gmx_fatal(FARGS,"Can not read Espresso files with particle property '%s'",word);
 484         }
 485         if (bFoundProp && p == espMOLECULE)
 486           bMol = TRUE;
 487         if (r == 3)
 488           level--;
 489       }
 490
 491       i = 0;
 492       while (level > 0 && (r=get_espresso_word(fp,word))) {
 493         if (r == 2) {
 494           level++;
 495         } else if (r == 3) {
 496           level--;
 497         }
 498         if (level == 2) {
 499           for(p=0; p<nprop; p++) {
 500             switch (prop[p]) {
 501             case espID:
 502               r = get_espresso_word(fp,word);
 503               /* Not used */
 504               break;
 505             case espPOS:
 506               for(m=0; m<3; m++) {
 507                 r = get_espresso_word(fp,word);
 508                 sscanf(word,"%lf",&d);
 509                 x[i][m] = d;
 510               }
 511               break;
 512             case espTYPE:
 513               r = get_espresso_word(fp,word);
 514               atoms->atom[i].type = strtol(word, NULL, 10);
 515               break;
 516             case espQ:
 517               r = get_espresso_word(fp,word);
 518               sscanf(word,"%lf",&d);
 519               atoms->atom[i].q = d;
 520               break;
 521             case espV:
 522               for(m=0; m<3; m++) {
 523                 r = get_espresso_word(fp,word);
 524                 sscanf(word,"%lf",&d);
 525                 v[i][m] = d;
 526               }
 527               break;
 528             case espF:
 529               for(m=0; m<3; m++) {
 530                 r = get_espresso_word(fp,word);
 531                 /* not used */
 532               }
 533               break;
 534             case espMOLECULE:
 535               r = get_espresso_word(fp,word);
 536               molnr = strtol(word, NULL, 10);
 537               if (i == 0 ||
 538                   atoms->resinfo[atoms->atom[i-1].resind].nr != molnr) {
 539                 atoms->atom[i].resind =
 540                   (i == 0 ? 0 : atoms->atom[i-1].resind+1);
 541                 atoms->resinfo[atoms->atom[i].resind].nr      = molnr;
 542                 atoms->resinfo[atoms->atom[i].resind].ic      = ' ';
 543                 atoms->resinfo[atoms->atom[i].resind].chainid = ' ';
 544         atoms->resinfo[atoms->atom[i].resind].chainnum = molnr; /* Not sure if this is right? */
 545               } else {
 546                 atoms->atom[i].resind = atoms->atom[i-1].resind;
 547               }
 548               break;
 549             }
 550           }
 551           /* Generate an atom name from the particle type */
 552           sprintf(buf,"T%d",atoms->atom[i].type);
 553           atoms->atomname[i] = put_symtab(symtab,buf);
 554           if (bMol) {
 555             if (i == 0 || atoms->atom[i].resind != atoms->atom[i-1].resind) {
 556               atoms->resinfo[atoms->atom[i].resind].name =
 557                 put_symtab(symtab,"MOL");
 558             }
 559           } else {
 560             /* Residue number is the atom number */
 561             atoms->atom[i].resind = i;
 562             /* Generate an residue name from the particle type */
 563             if (atoms->atom[i].type < 26) {
 564               sprintf(buf,"T%c",'A'+atoms->atom[i].type);
 565             } else {
 566               sprintf(buf,"T%c%c",
 567                       'A'+atoms->atom[i].type/26,'A'+atoms->atom[i].type%26);
 568             }
 569             t_atoms_set_resinfo(atoms,i,symtab,buf,i,' ',0,' ');
 570           }
 571
 572           if (r == 3)
 573             level--;
 574           i++;
 575         }
 576       }
 577       atoms->nres = atoms->nr;
 578
 579       if (i != atoms->nr) {
 580         gmx_fatal(FARGS,"Internal inconsistency in Espresso routines, read %d atoms, expected %d atoms",i,atoms->nr);
 581       }
 582     } else if (level==1 && strcmp(word,"variable")==0 && !bFoundVariable) {
 583       bFoundVariable = TRUE;
 584       level += check_open_parenthesis(fp,r,infile,"variable");
 585       while (level==2 && (r=get_espresso_word(fp,word))) {
 586         if (level==2 && strcmp(word,"box_l") == 0) {
 587           for(m=0; m<3; m++) {
 588             r = get_espresso_word(fp,word);
 589             sscanf(word,"%lf",&d);
 590             box[m][m] = d;
 591           }
 592           level += check_close_parenthesis(fp,r,infile,"box_l");
 593         }
 594       }
 595     } else if (r == 2) {
 596       level++;
 597     } else if (r == 3) {
 598       level--;
 599     }
 600   }
 601
 602   if (!bFoundParticles) {
 603     fprintf(stderr,"Did not find a particles section in Espresso file '%s'\n",
 604             infile);
 605   }
 606
 607   gmx_fio_fclose(fp);
 608 }
 609
 610 static int get_espresso_coordnum(const char *infile)
 611 {
 612   FILE *fp;
 613   char word[STRLEN];
 614   int  natoms,level,r;
 615   gmx_bool bFoundParticles;
 616
 617   natoms = 0;
 618
 619   fp = gmx_fio_fopen(infile,"r");
 620
 621   bFoundParticles = FALSE;
 622   level = 0;
 623   while ((r=get_espresso_word(fp,word)) && !bFoundParticles) {
 624     if (level==1 && strcmp(word,"particles")==0 && !bFoundParticles) {
 625       bFoundParticles = TRUE;
 626       level += check_open_parenthesis(fp,r,infile,"particles");
 627       while (level > 0 && (r=get_espresso_word(fp,word))) {
 628         if (r == 2) {
 629           level++;
 630           if (level == 2)
 631             natoms++;
 632         } else if (r == 3) {
 633           level--;
 634         }
 635       }
 636     } else if (r == 2) {
 637       level++;
 638     } else if (r == 3) {
 639       level--;
 640     }
 641   }
 642   if (!bFoundParticles) {
 643     fprintf(stderr,"Did not find a particles section in Espresso file '%s'\n",
 644             infile);
 645   }
 646
 647   gmx_fio_fclose(fp);
 648
 649   return natoms;
 650 }
 651
 652 static void write_espresso_conf_indexed(FILE *out,const char *title,
 653                                         t_atoms *atoms,int nx,atom_id *index,
 654                                         rvec *x,rvec *v,matrix box)
 655 {
 656   int i,j;
 657
 658   fprintf(out,"# %s\n",title);
 659   if (TRICLINIC(box)) {
 660     gmx_warning("The Espresso format does not support triclinic unit-cells");
 661   }
 662   fprintf(out,"{variable {box_l %f %f %f}}\n",box[0][0],box[1][1],box[2][2]);
 663
 664   fprintf(out,"{particles {id pos type q%s}\n",v ? " v" : "");
 665   for(i=0; i<nx; i++) {
 666     if (index)
 667       j = index[i];
 668     else
 669       j = i;
 670     fprintf(out,"\t{%d %f %f %f %d %g",
 671             j,x[j][XX],x[j][YY],x[j][ZZ],
 672             atoms->atom[j].type,atoms->atom[j].q);
 673     if (v)
 674       fprintf(out," %f %f %f",v[j][XX],v[j][YY],v[j][ZZ]);
 675     fprintf(out,"}\n");
 676   }
 677   fprintf(out,"}\n");
 678 }
 679
 680 static void get_coordnum_fp (FILE *in, char *title, int *natoms)
 681 {
 682   char line[STRLEN+1];
 683
 684   fgets2 (title,STRLEN,in);
 685   fgets2 (line,STRLEN,in);
 686   if (sscanf (line,"%d",natoms) != 1) {
 687     gmx_fatal(FARGS,"gro file does not have the number of atoms on the second line");
 688   }
 689 }
 690
 691 static void get_coordnum (const char *infile,int *natoms)
 692 {
 693   FILE *in;
 694   char title[STRLEN];
 695
 696   in=gmx_fio_fopen(infile,"r");
 697   get_coordnum_fp(in,title,natoms);
 698   gmx_fio_fclose (in);
 699 }
 700
 701 static gmx_bool get_w_conf(FILE *in, const char *infile, char *title,
 702                            t_symtab *symtab, t_atoms *atoms, int *ndec,
 703                            rvec x[], rvec *v, matrix box)
 704 {
 705   char   name[6];
 706   char   line[STRLEN+1],*ptr;
 707   char   buf[256];
 708   double x1,y1,z1,x2,y2,z2;
 709   rvec   xmin,xmax;
 710   int    natoms,i,m,resnr,newres,oldres,ddist,c;
 711   gmx_bool   bFirst,bVel;
 712   char   *p1,*p2,*p3;
 713
 714   newres  = -1;
 715   oldres  = NOTSET; /* Unlikely number for the first residue! */
 716   ddist   = 0;
 717
 718   /* Read the title and number of atoms */
 719   get_coordnum_fp(in,title,&natoms);
 720
 721   if (natoms > atoms->nr)
 722     gmx_fatal(FARGS,"gro file contains more atoms (%d) than expected (%d)",
 723                 natoms,atoms->nr);
 724   else if (natoms <  atoms->nr)
 725     fprintf(stderr,"Warning: gro file contains less atoms (%d) than expected"
 726             " (%d)\n",natoms,atoms->nr);
 727
 728   bFirst=TRUE;
 729
 730   bVel = FALSE;
 731
 732   /* just pray the arrays are big enough */
 733   for (i=0; (i < natoms) ; i++) {
 734     if ((fgets2 (line,STRLEN,in)) == NULL) {
 735       unexpected_eof(infile,i+2);
 736     }
 737     if (strlen(line) < 39)
 738       gmx_fatal(FARGS,"Invalid line in %s for atom %d:\n%s",infile,i+1,line);
 739
 740     /* determine read precision from distance between periods
 741        (decimal points) */
 742     if (bFirst) {
 743       bFirst=FALSE;
 744       p1=strchr(line,'.');
 745       if (p1 == NULL)
 746         gmx_fatal(FARGS,"A coordinate in file %s does not contain a '.'",infile);
 747       p2=strchr(&p1[1],'.');
 748       if (p2 == NULL)
 749         gmx_fatal(FARGS,"A coordinate in file %s does not contain a '.'",infile);
 750       ddist = p2 - p1;
 751       *ndec = ddist - 5;
 752
 753       p3=strchr(&p2[1],'.');
 754       if (p3 == NULL)
 755         gmx_fatal(FARGS,"A coordinate in file %s does not contain a '.'",infile);
 756
 757       if (p3 - p2 != ddist)
 758         gmx_fatal(FARGS,"The spacing of the decimal points in file %s is not consistent for x, y and z",infile);
 759     }
 760
 761     /* residue number*/
 762     memcpy(name,line,5);
 763     name[5]='\0';
 764     sscanf(name,"%d",&resnr);
 765     memcpy(name,line+5,5);
 766     name[5]='\0';
 767     if (resnr != oldres) {
 768       oldres = resnr;
 769       newres++;
 770       if (newres >= natoms)
 771         gmx_fatal(FARGS,"More residues than atoms in %s (natoms = %d)",
 772                     infile,natoms);
 773       atoms->atom[i].resind = newres;
 774       t_atoms_set_resinfo(atoms,i,symtab,name,resnr,' ',0,' ');
 775     } else {
 776       atoms->atom[i].resind = newres;
 777     }
 778
 779     /* atomname */
 780     memcpy(name,line+10,5);
 781     atoms->atomname[i]=put_symtab(symtab,name);
 782
 783     /* eventueel controle atomnumber met i+1 */
 784
 785     /* coordinates (start after residue shit) */
 786     ptr = line + 20;
 787     /* Read fixed format */
 788     for(m=0; m<DIM; m++) {
 789       for(c=0; (c<ddist && ptr[0]); c++) {
 790         buf[c] = ptr[0];
 791         ptr++;
 792       }
 793       buf[c] = '\0';
 794       if (sscanf (buf,"%lf %lf",&x1,&x2) != 1) {
 795         gmx_fatal(FARGS,"Something is wrong in the coordinate formatting of file %s. Note that gro is fixed format (see the manual)",infile);
 796       } else {
 797         x[i][m] = x1;
 798       }
 799     }
 800
 801     /* velocities (start after residues and coordinates) */
 802     if (v) {
 803       /* Read fixed format */
 804       for(m=0; m<DIM; m++) {
 805         for(c=0; (c<ddist && ptr[0]); c++) {
 806           buf[c] = ptr[0];
 807           ptr++;
 808         }
 809         buf[c] = '\0';
 810         if (sscanf (buf,"%lf",&x1) != 1) {
 811           v[i][m] = 0;
 812         } else {
 813           v[i][m] = x1;
 814           bVel = TRUE;
 815         }
 816       }
 817     }
 818   }
 819   atoms->nres = newres + 1;
 820
 821   /* box */
 822   fgets2 (line,STRLEN,in);
 823   if (sscanf (line,"%lf%lf%lf",&x1,&y1,&z1) != 3) {
 824     gmx_warning("Bad box in file %s",infile);
 825
 826     /* Generate a cubic box */
 827     for(m=0; (m<DIM); m++)
 828       xmin[m]=xmax[m]=x[0][m];
 829     for(i=1; (i<atoms->nr); i++)
 830       for(m=0; (m<DIM); m++) {
 831         xmin[m]=min(xmin[m],x[i][m]);
 832         xmax[m]=max(xmax[m],x[i][m]);
 833       }
 834     for (i=0; i<DIM; i++) for (m=0; m<DIM; m++) box[i][m]=0.0;
 835     for(m=0; (m<DIM); m++)
 836       box[m][m]=(xmax[m]-xmin[m]);
 837     fprintf(stderr,"Generated a cubic box %8.3f x %8.3f x %8.3f\n",
 838             box[XX][XX],box[YY][YY],box[ZZ][ZZ]);
 839   }
 840   else {
 841     /* We found the first three values, the diagonal elements */
 842     box[XX][XX]=x1;
 843     box[YY][YY]=y1;
 844     box[ZZ][ZZ]=z1;
 845     if (sscanf (line,"%*f%*f%*f%lf%lf%lf%lf%lf%lf",
 846                 &x1,&y1,&z1,&x2,&y2,&z2) != 6)
 847       x1=y1=z1=x2=y2=z2=0.0;
 848     box[XX][YY] = x1;
 849     box[XX][ZZ] = y1;
 850     box[YY][XX] = z1;
 851     box[YY][ZZ] = x2;
 852     box[ZZ][XX] = y2;
 853     box[ZZ][YY] = z2;
 854   }
 855
 856   return bVel;
 857 }
 858
 859 static void read_whole_conf(const char *infile,char *title,
 860                             t_atoms *atoms, rvec x[],rvec *v, matrix box)
 861 {
 862   FILE   *in;
 863   int    ndec;
 864   t_symtab symtab;
 865
 866   /* open file */
 867   in=gmx_fio_fopen(infile,"r");
 868
 869   open_symtab(&symtab);
 870   get_w_conf(in, infile, title, &symtab, atoms, &ndec, x, v, box);
 871   /* We can't free the symbols, as they are still used in atoms, so
 872    * the only choice is to leak them. */
 873   free_symtab(&symtab);
 874
 875   gmx_fio_fclose(in);
 876 }
 877
 878 gmx_bool gro_next_x_or_v(FILE *status,t_trxframe *fr)
 879 {
 880   t_atoms atoms;
 881   t_symtab symtab;
 882   char    title[STRLEN],*p;
 883   double  tt;
 884   int     ndec=0,i;
 885
 886   if (gmx_eof(status))
 887     return FALSE;
 888
 889   open_symtab(&symtab);
 890   atoms.nr=fr->natoms;
 891   snew(atoms.atom,fr->natoms);
 892   atoms.nres=fr->natoms;
 893   snew(atoms.resinfo,fr->natoms);
 894   snew(atoms.atomname,fr->natoms);
 895
 896   fr->bV = get_w_conf(status,title,title,&symtab,&atoms,&ndec,fr->x,fr->v,fr->box);
 897   fr->bPrec = TRUE;
 898   fr->prec = 1;
 899   /* prec = 10^ndec: */
 900   for(i=0; i<ndec; i++)
 901     fr->prec *= 10;
 902   fr->title = title;
 903   fr->bTitle = TRUE;
 904   fr->bX = TRUE;
 905   fr->bBox = TRUE;
 906
 907   sfree(atoms.atom);
 908   sfree(atoms.resinfo);
 909   sfree(atoms.atomname);
 910   done_symtab(&symtab);
 911
 912   if ((p=strstr(title,"t=")) != NULL) {
 913     p+=2;
 914     if (sscanf(p,"%lf",&tt)==1) {
 915       fr->time = tt;
 916       fr->bTime = TRUE;
 917     } else {
 918       fr->time = 0;
 919       fr->bTime = FALSE;
 920     }
 921   }
 922
 923   if (atoms.nr != fr->natoms)
 924     gmx_fatal(FARGS,"Number of atoms in gro frame (%d) doesn't match the number in the previous frame (%d)",atoms.nr,fr->natoms);
 925
 926   return TRUE;
 927 }
 928
 929 int gro_first_x_or_v(FILE *status,t_trxframe *fr)
 930 {
 931   char title[STRLEN];
 932
 933   frewind(status);
 934   fprintf(stderr,"Reading frames from gro file");
 935   get_coordnum_fp(status, title, &fr->natoms);
 936   frewind(status);
 937   fprintf(stderr," '%s', %d atoms.\n",title, fr->natoms);
 938   fr->bTitle = TRUE;
 939   fr->title = title;
 940   if (fr->natoms==0)
 941     gmx_file("No coordinates in gro file");
 942
 943   snew(fr->x,fr->natoms);
 944   snew(fr->v,fr->natoms);
 945   gro_next_x_or_v(status, fr);
 946
 947   return fr->natoms;
 948 }
 949
 950 static void make_hconf_format(int pr,gmx_bool bVel,char format[])
 951 {
 952   int l,vpr;
 953
 954   /* build format string for printing,
 955      something like "%8.3f" for x and "%8.4f" for v */
 956   if (pr<0)
 957     pr=0;
 958   if (pr>30)
 959     pr=30;
 960   l=pr+5;
 961   vpr=pr+1;
 962   if (bVel)
 963     sprintf(format,"%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
 964             l,pr,l,pr,l,pr,l,vpr,l,vpr,l,vpr);
 965   else
 966     sprintf(format,"%%%d.%df%%%d.%df%%%d.%df\n",l,pr,l,pr,l,pr);
 967
 968 }
 969
 970 static void write_hconf_box(FILE *out,int pr,matrix box)
 971 {
 972   char format[100];
 973   int l;
 974
 975   if (pr<5)
 976     pr=5;
 977   l=pr+5;
 978
 979   if (box[XX][YY] || box[XX][ZZ] || box[YY][XX] || box[YY][ZZ] ||
 980       box[ZZ][XX] || box[ZZ][YY]) {
 981     sprintf(format,"%%%d.%df%%%d.%df%%%d.%df"
 982             "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
 983             l,pr,l,pr,l,pr,l,pr,l,pr,l,pr,l,pr,l,pr,l,pr);
 984     fprintf(out,format,
 985             box[XX][XX],box[YY][YY],box[ZZ][ZZ],
 986             box[XX][YY],box[XX][ZZ],box[YY][XX],
 987             box[YY][ZZ],box[ZZ][XX],box[ZZ][YY]);
 988   } else {
 989     sprintf(format,"%%%d.%df%%%d.%df%%%d.%df\n",l,pr,l,pr,l,pr);
 990     fprintf(out,format,
 991             box[XX][XX],box[YY][YY],box[ZZ][ZZ]);
 992   }
 993 }
 994
 995 void write_hconf_indexed_p(FILE *out,const char *title,t_atoms *atoms,
 996                            int nx,const atom_id index[], int pr,
 997                            rvec *x,rvec *v,matrix box)
 998 {
 999   char resnm[6],nm[6],format[100];
1000   int  ai,i,resind,resnr;
1001
1002   bromacs(format,99);
1003   fprintf (out,"%s\n",(title && title[0])?title:format);
1004   fprintf (out,"%5d\n",nx);
1005
1006   make_hconf_format(pr,v!=NULL,format);
1007
1008   for (i=0; (i<nx); i++) {
1009     ai=index[i];
1010
1011     resind = atoms->atom[ai].resind;
1012     strncpy(resnm," ??? ",sizeof(resnm)-1);
1013     if (resind < atoms->nres) {
1014       strncpy(resnm,*atoms->resinfo[resind].name,sizeof(resnm)-1);
1015       resnr = atoms->resinfo[resind].nr;
1016     } else {
1017       strncpy(resnm," ??? ",sizeof(resnm)-1);
1018       resnr = resind + 1;
1019     }
1020
1021     if (atoms->atom)
1022       strncpy(nm,*atoms->atomname[ai],sizeof(nm)-1);
1023     else
1024       strncpy(nm," ??? ",sizeof(nm)-1);
1025
1026     fprintf(out,"%5d%-5.5s%5.5s%5d",resnr%100000,resnm,nm,(ai+1)%100000);
1027     /* next fprintf uses built format string */
1028     if (v)
1029       fprintf(out,format,
1030               x[ai][XX], x[ai][YY], x[ai][ZZ], v[ai][XX],v[ai][YY],v[ai][ZZ]);
1031     else
1032       fprintf(out,format,
1033               x[ai][XX], x[ai][YY], x[ai][ZZ]);
1034   }
1035
1036   write_hconf_box(out,pr,box);
1037
1038   fflush(out);
1039 }
1040
1041 static void write_hconf_mtop(FILE *out,const char *title,gmx_mtop_t *mtop,
1042                              int pr,
1043                              rvec *x,rvec *v,matrix box)
1044 {
1045   char format[100];
1046   int  i,resnr;
1047   gmx_mtop_atomloop_all_t aloop;
1048   t_atom *atom;
1049   char *atomname,*resname;
1050
1051   bromacs(format,99);
1052   fprintf (out,"%s\n",(title && title[0])?title:format);
1053   fprintf (out,"%5d\n",mtop->natoms);
1054
1055   make_hconf_format(pr,v!=NULL,format);
1056
1057   aloop = gmx_mtop_atomloop_all_init(mtop);
1058   while (gmx_mtop_atomloop_all_next(aloop,&i,&atom)) {
1059     gmx_mtop_atomloop_all_names(aloop,&atomname,&resnr,&resname);
1060
1061     fprintf(out,"%5d%-5.5s%5.5s%5d",
1062             resnr%100000,resname,atomname,(i+1)%100000);
1063     /* next fprintf uses built format string */
1064     if (v)
1065       fprintf(out,format,
1066               x[i][XX], x[i][YY], x[i][ZZ], v[i][XX],v[i][YY],v[i][ZZ]);
1067     else
1068       fprintf(out,format,
1069               x[i][XX], x[i][YY], x[i][ZZ]);
1070   }
1071
1072   write_hconf_box(out,pr,box);
1073
1074   fflush(out);
1075 }
1076
1077 void write_hconf_p(FILE *out,const char *title,t_atoms *atoms, int pr,
1078                    rvec *x,rvec *v,matrix box)
1079 {
1080   atom_id *aa;
1081   int     i;
1082
1083   snew(aa,atoms->nr);
1084   for(i=0; (i<atoms->nr); i++)
1085     aa[i]=i;
1086   write_hconf_indexed_p(out,title,atoms,atoms->nr,aa,pr,x,v,box);
1087   sfree(aa);
1088 }
1089
1090 void write_conf_p(const char *outfile, const char *title,
1091                   t_atoms *atoms, int pr,
1092                   rvec *x, rvec *v,matrix box)
1093 {
1094   FILE *out;
1095
1096   out=gmx_fio_fopen(outfile,"w");
1097   write_hconf_p(out,title,atoms,pr,x,v,box);
1098
1099   gmx_fio_fclose (out);
1100 }
1101
1102 static void write_conf(const char *outfile, const char *title, t_atoms *atoms,
1103                        rvec *x, rvec *v,matrix box)
1104 {
1105   write_conf_p(outfile, title, atoms, 3, x, v, box);
1106 }
1107
1108 void write_sto_conf_indexed(const char *outfile,const char *title,
1109                             t_atoms *atoms,
1110                             rvec x[],rvec *v,int ePBC,matrix box,
1111                             atom_id nindex,atom_id index[])
1112 {
1113   FILE       *out;
1114   int        ftp;
1115   t_trxframe fr;
1116
1117   ftp=fn2ftp(outfile);
1118   switch (ftp) {
1119   case efGRO:
1120     out=gmx_fio_fopen(outfile,"w");
1121     write_hconf_indexed_p(out, title, atoms, nindex, index, 3, x, v, box);
1122     gmx_fio_fclose(out);
1123     break;
1124   case efG96:
1125     clear_trxframe(&fr,TRUE);
1126     fr.bTitle = TRUE;
1127     fr.title = title;
1128     fr.natoms = atoms->nr;
1129     fr.bAtoms = TRUE;
1130     fr.atoms = atoms;
1131     fr.bX = TRUE;
1132     fr.x = x;
1133     if (v) {
1134       fr.bV = TRUE;
1135       fr.v = v;
1136     }
1137     fr.bBox = TRUE;
1138     copy_mat(box,fr.box);
1139     out=gmx_fio_fopen(outfile,"w");
1140     write_g96_conf(out, &fr, nindex, index);
1141     gmx_fio_fclose(out);
1142     break;
1143   case efPDB:
1144   case efBRK:
1145   case efENT:
1146   case efPQR:
1147     out=gmx_fio_fopen(outfile,"w");
1148     write_pdbfile_indexed(out,title,atoms,x,ePBC,box,' ',-1,nindex,index,NULL,TRUE);
1149     gmx_fio_fclose(out);
1150     break;
1151   case efESP:
1152     out=gmx_fio_fopen(outfile,"w");
1153     write_espresso_conf_indexed(out, title, atoms, nindex, index, x, v, box);
1154     gmx_fio_fclose(out);
1155     break;
1156   case efTPR:
1157   case efTPB:
1158   case efTPA:
1159     gmx_fatal(FARGS,"Sorry, can not write a topology to %s",outfile);
1160     break;
1161   default:
1162     gmx_incons("Not supported in write_sto_conf_indexed");
1163   }
1164 }
1165
1166 static void write_xyz_conf(const char *outfile,const char *title,
1167                            t_atoms *atoms,rvec *x)
1168 {
1169   FILE *fp;
1170   int i,anr;
1171   real value;
1172   char *ptr,*name;
1173   gmx_atomprop_t aps=gmx_atomprop_init();
1174
1175   fp = gmx_fio_fopen(outfile,"w");
1176   fprintf(fp,"%3d\n",atoms->nr);
1177   fprintf(fp,"%s\n",title);
1178   for(i=0; (i<atoms->nr); i++) {
1179     anr  = atoms->atom[i].atomnumber;
1180     name = *atoms->atomname[i];
1181     if (anr == NOTSET) {
1182       if (gmx_atomprop_query(aps,epropElement,"???",name,&value))
1183         anr = gmx_nint(value);
1184     }
1185     if ((ptr = gmx_atomprop_element(aps,anr)) == NULL)
1186       ptr = name;
1187     fprintf(fp,"%3s%10.5f%10.5f%10.5f\n",ptr,
1188             10*x[i][XX],10*x[i][YY],10*x[i][ZZ]);
1189   }
1190   gmx_fio_fclose(fp);
1191   gmx_atomprop_destroy(aps);
1192 }
1193
1194 void write_sto_conf(const char *outfile,const char *title,t_atoms *atoms,
1195                     rvec x[],rvec *v,int ePBC,matrix box)
1196 {
1197   FILE       *out;
1198   int        ftp;
1199   t_trxframe fr;
1200
1201   ftp=fn2ftp(outfile);
1202   switch (ftp) {
1203   case efGRO:
1204     write_conf(outfile, title, atoms, x, v, box);
1205     break;
1206   case efG96:
1207     clear_trxframe(&fr,TRUE);
1208     fr.bTitle = TRUE;
1209     fr.title = title;
1210     fr.natoms = atoms->nr;
1211     fr.bAtoms = TRUE;
1212     fr.atoms = atoms;
1213     fr.bX = TRUE;
1214     fr.x = x;
1215     if (v) {
1216       fr.bV = TRUE;
1217       fr.v = v;
1218     }
1219     fr.bBox = TRUE;
1220     copy_mat(box,fr.box);
1221     out=gmx_fio_fopen(outfile,"w");
1222     write_g96_conf(out, &fr, -1, NULL);
1223     gmx_fio_fclose(out);
1224     break;
1225   case efXYZ:
1226     write_xyz_conf(outfile,(strlen(title) > 0) ? title : outfile,atoms,x);
1227     break;
1228   case efPDB:
1229   case efBRK:
1230   case efENT:
1231     out=gmx_fio_fopen(outfile,"w");
1232     write_pdbfile(out, title, atoms, x, ePBC, box, ' ', -1,NULL,TRUE);
1233     gmx_fio_fclose(out);
1234     break;
1235   case efESP:
1236     out=gmx_fio_fopen(outfile,"w");
1237     write_espresso_conf_indexed(out, title, atoms, atoms->nr, NULL, x, v, box);
1238     gmx_fio_fclose(out);
1239     break;
1240   case efTPR:
1241   case efTPB:
1242   case efTPA:
1243     gmx_fatal(FARGS,"Sorry, can not write a topology to %s",outfile);
1244     break;
1245   default:
1246     gmx_incons("Not supported in write_sto_conf");
1247   }
1248 }
1249
1250 void write_sto_conf_mtop(const char *outfile,const char *title,
1251                          gmx_mtop_t *mtop,
1252                          rvec x[],rvec *v,int ePBC,matrix box)
1253 {
1254   int  ftp;
1255   FILE *out;
1256   t_atoms atoms;
1257
1258   ftp=fn2ftp(outfile);
1259   switch (ftp) {
1260   case efGRO:
1261     out = gmx_fio_fopen(outfile,"w");
1262     write_hconf_mtop(out,title,mtop,3,x,v,box);
1263     gmx_fio_fclose(out);
1264     break;
1265   default:
1266     /* This is a brute force approach which requires a lot of memory.
1267      * We should implement mtop versions of all writing routines.
1268      */
1269     atoms = gmx_mtop_global_atoms(mtop);
1270
1271     write_sto_conf(outfile,title,&atoms,x,v,ePBC,box);
1272
1273     done_atom(&atoms);
1274     break;
1275   }
1276 }
1277
1278 static int get_xyz_coordnum(const char *infile)
1279 {
1280   FILE *fp;
1281   int n;
1282
1283   fp = gmx_fio_fopen(infile,"r");
1284   if (fscanf(fp,"%d",&n) != 1)
1285     gmx_fatal(FARGS,"Can not read number of atoms from %s",infile);
1286   gmx_fio_fclose(fp);
1287
1288   return n;
1289 }
1290
1291 static void read_xyz_conf(const char *infile,char *title,
1292                           t_atoms *atoms,rvec *x)
1293 {
1294   FILE   *fp;
1295   int    i,n;
1296   double xx,yy,zz;
1297   t_symtab *tab;
1298   char atomnm[32],buf[STRLEN];
1299
1300   snew(tab,1);
1301   fp = gmx_fio_fopen(infile,"r");
1302   fgets2(buf,STRLEN-1,fp);
1303   if (sscanf(buf,"%d",&n) != 1)
1304     gmx_fatal(FARGS,"Can not read number of atoms from %s",infile);
1305   fgets2(buf,STRLEN-1,fp);
1306   strcpy(title,buf);
1307   for(i=0; (i<n); i++) {
1308     fgets2(buf,STRLEN-1,fp);
1309     if (sscanf(buf,"%s%lf%lf%lf",atomnm,&xx,&yy,&zz) != 4)
1310       gmx_fatal(FARGS,"Can not read coordinates from %s",infile);
1311     atoms->atomname[i] = put_symtab(tab,atomnm);
1312     x[i][XX] = xx*0.1;
1313     x[i][YY] = yy*0.1;
1314     x[i][ZZ] = zz*0.1;
1315   }
1316   gmx_fio_fclose(fp);
1317 }
1318
1319 void get_stx_coordnum(const char *infile,int *natoms)
1320 {
1321   FILE *in;
1322   int ftp,tpxver,tpxgen;
1323   t_trxframe fr;
1324   char g96_line[STRLEN+1];
1325
1326   ftp=fn2ftp(infile);
1327   range_check(ftp,0,efNR);
1328   switch (ftp) {
1329   case efGRO:
1330     get_coordnum(infile, natoms);
1331     break;
1332   case efG96:
1333     in=gmx_fio_fopen(infile,"r");
1334     fr.title = NULL;
1335     fr.natoms = -1;
1336     fr.atoms = NULL;
1337     fr.x = NULL;
1338     fr.v = NULL;
1339     fr.f = NULL;
1340     *natoms=read_g96_conf(in,infile,&fr,g96_line);
1341     gmx_fio_fclose(in);
1342     break;
1343   case efXYZ:
1344     *natoms = get_xyz_coordnum(infile);
1345     break;
1346   case efPDB:
1347   case efBRK:
1348   case efENT:
1349     in=gmx_fio_fopen(infile,"r");
1350     get_pdb_coordnum(in, natoms);
1351     gmx_fio_fclose(in);
1352     break;
1353   case efESP:
1354     *natoms = get_espresso_coordnum(infile);
1355     break;
1356   case efTPA:
1357   case efTPB:
1358   case efTPR: {
1359     t_tpxheader tpx;
1360
1361     read_tpxheader(infile,&tpx,TRUE,&tpxver,&tpxgen);
1362     *natoms = tpx.natoms;
1363     break;
1364   }
1365   default:
1366     gmx_fatal(FARGS,"File type %s not supported in get_stx_coordnum",
1367               ftp2ext(ftp));
1368   }
1369 }
1370
1371 void read_stx_conf(const char *infile,char *title,t_atoms *atoms,
1372                    rvec x[],rvec *v,int *ePBC,matrix box)
1373 {
1374   FILE       *in;
1375   char       buf[256];
1376   gmx_mtop_t *mtop;
1377   t_topology top;
1378   t_trxframe fr;
1379   int        i,ftp,natoms;
1380   real       d;
1381   char       g96_line[STRLEN+1];
1382
1383   if (atoms->nr == 0)
1384     fprintf(stderr,"Warning: Number of atoms in %s is 0\n",infile);
1385   else if (atoms->atom == NULL)
1386     gmx_mem("Uninitialized array atom");
1387
1388   if (ePBC)
1389     *ePBC = -1;
1390
1391   ftp=fn2ftp(infile);
1392   switch (ftp) {
1393   case efGRO:
1394     read_whole_conf(infile, title, atoms, x, v, box);
1395     break;
1396   case efXYZ:
1397     read_xyz_conf(infile,title,atoms,x);
1398     break;
1399   case efG96:
1400     fr.title  = NULL;
1401     fr.natoms = atoms->nr;
1402     fr.atoms = atoms;
1403     fr.x = x;
1404     fr.v = v;
1405     fr.f = NULL;
1406     in = gmx_fio_fopen(infile,"r");
1407     read_g96_conf(in, infile, &fr, g96_line);
1408     gmx_fio_fclose(in);
1409     copy_mat(fr.box,box);
1410     strncpy(title, fr.title, STRLEN);
1411     break;
1412   case efPDB:
1413   case efBRK:
1414   case efENT:
1415     read_pdb_conf(infile, title, atoms, x, ePBC, box, TRUE, NULL);
1416     break;
1417   case efESP:
1418     read_espresso_conf(infile,atoms,x,v,box);
1419     break;
1420   case efTPR:
1421   case efTPB:
1422   case efTPA:
1423     snew(mtop,1);
1424     i = read_tpx(infile,NULL,box,&natoms,x,v,NULL,mtop);
1425     if (ePBC)
1426       *ePBC = i;
1427
1428     strcpy(title,*(mtop->name));
1429
1430     /* Free possibly allocated memory */
1431     done_atom(atoms);
1432
1433     *atoms = gmx_mtop_global_atoms(mtop);
1434     top = gmx_mtop_t_to_t_topology(mtop);
1435     tpx_make_chain_identifiers(atoms,&top.mols);
1436
1437     sfree(mtop);
1438     /* The strings in the symtab are still in use in the returned t_atoms
1439      * structure, so we should not free them. But there is no place to put the
1440      * symbols; the only choice is to leak the memory...
1441      * So we clear the symbol table before freeing the topology structure. */
1442     free_symtab(&top.symtab);
1443     done_top(&top);
1444
1445     break;
1446   default:
1447     gmx_incons("Not supported in read_stx_conf");
1448   }
1449 }
1450