src/gromacs/fileio/confio.c

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