src/gromacs/fileio/confio.cpp

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