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