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