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