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