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, char *title,
 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     // TODO: The code does not understand titles it writes...
 619     title[0] = '\0';
 620
 621     clear_mat(box);
 622
 623     fp = gmx_fio_fopen(infile, "r");
 624
 625     bFoundParticles = FALSE;
 626     bFoundVariable  = FALSE;
 627     bMol            = FALSE;
 628     level           = 0;
 629     while ((r = get_espresso_word(fp, word)))
 630     {
 631         if (level == 1 && std::strcmp(word, "particles") == 0 && !bFoundParticles)
 632         {
 633             bFoundParticles = TRUE;
 634             level          += check_open_parenthesis(fp, r, infile, "particles");
 635             nprop           = 0;
 636             while (level == 2 && (r = get_espresso_word(fp, word)))
 637             {
 638                 bFoundProp = FALSE;
 639                 for (p = 0; p < espNR; p++)
 640                 {
 641                     if (strcmp(word, esp_prop[p]) == 0)
 642                     {
 643                         bFoundProp    = TRUE;
 644                         prop[nprop++] = p;
 645                         /* printf("  prop[%d] = %s\n",nprop-1,esp_prop[prop[nprop-1]]); */
 646                     }
 647                 }
 648                 if (!bFoundProp && word[0] != '}')
 649                 {
 650                     gmx_fatal(FARGS, "Can not read Espresso files with particle property '%s'", word);
 651                 }
 652                 if (bFoundProp && p == espMOLECULE)
 653                 {
 654                     bMol = TRUE;
 655                 }
 656                 if (r == 3)
 657                 {
 658                     level--;
 659                 }
 660             }
 661
 662             i = 0;
 663             while (level > 0 && (r = get_espresso_word(fp, word)))
 664             {
 665                 if (r == 2)
 666                 {
 667                     level++;
 668                 }
 669                 else if (r == 3)
 670                 {
 671                     level--;
 672                 }
 673                 if (level == 2)
 674                 {
 675                     for (p = 0; p < nprop; p++)
 676                     {
 677                         switch (prop[p])
 678                         {
 679                             case espID:
 680                                 r = get_espresso_word(fp, word);
 681                                 /* Not used */
 682                                 break;
 683                             case espPOS:
 684                                 for (m = 0; m < 3; m++)
 685                                 {
 686                                     r = get_espresso_word(fp, word);
 687                                     sscanf(word, "%lf", &d);
 688                                     x[i][m] = d;
 689                                 }
 690                                 break;
 691                             case espTYPE:
 692                                 r                   = get_espresso_word(fp, word);
 693                                 atoms->atom[i].type = std::strtol(word, NULL, 10);
 694                                 break;
 695                             case espQ:
 696                                 r = get_espresso_word(fp, word);
 697                                 sscanf(word, "%lf", &d);
 698                                 atoms->atom[i].q = d;
 699                                 break;
 700                             case espV:
 701                                 for (m = 0; m < 3; m++)
 702                                 {
 703                                     r = get_espresso_word(fp, word);
 704                                     sscanf(word, "%lf", &d);
 705                                     v[i][m] = d;
 706                                 }
 707                                 break;
 708                             case espF:
 709                                 for (m = 0; m < 3; m++)
 710                                 {
 711                                     r = get_espresso_word(fp, word);
 712                                     /* not used */
 713                                 }
 714                                 break;
 715                             case espMOLECULE:
 716                                 r     = get_espresso_word(fp, word);
 717                                 molnr = std::strtol(word, NULL, 10);
 718                                 if (i == 0 ||
 719                                     atoms->resinfo[atoms->atom[i-1].resind].nr != molnr)
 720                                 {
 721                                     atoms->atom[i].resind =
 722                                         (i == 0 ? 0 : atoms->atom[i-1].resind+1);
 723                                     atoms->resinfo[atoms->atom[i].resind].nr       = molnr;
 724                                     atoms->resinfo[atoms->atom[i].resind].ic       = ' ';
 725                                     atoms->resinfo[atoms->atom[i].resind].chainid  = ' ';
 726                                     atoms->resinfo[atoms->atom[i].resind].chainnum = molnr; /* Not sure if this is right? */
 727                                 }
 728                                 else
 729                                 {
 730                                     atoms->atom[i].resind = atoms->atom[i-1].resind;
 731                                 }
 732                                 break;
 733                         }
 734                     }
 735                     /* Generate an atom name from the particle type */
 736                     sprintf(buf, "T%d", atoms->atom[i].type);
 737                     atoms->atomname[i] = put_symtab(symtab, buf);
 738                     if (bMol)
 739                     {
 740                         if (i == 0 || atoms->atom[i].resind != atoms->atom[i-1].resind)
 741                         {
 742                             atoms->resinfo[atoms->atom[i].resind].name =
 743                                 put_symtab(symtab, "MOL");
 744                         }
 745                     }
 746                     else
 747                     {
 748                         /* Residue number is the atom number */
 749                         atoms->atom[i].resind = i;
 750                         /* Generate an residue name from the particle type */
 751                         if (atoms->atom[i].type < 26)
 752                         {
 753                             sprintf(buf, "T%c", 'A'+atoms->atom[i].type);
 754                         }
 755                         else
 756                         {
 757                             sprintf(buf, "T%c%c",
 758                                     'A'+atoms->atom[i].type/26, 'A'+atoms->atom[i].type%26);
 759                         }
 760                         t_atoms_set_resinfo(atoms, i, symtab, buf, i, ' ', 0, ' ');
 761                     }
 762
 763                     if (r == 3)
 764                     {
 765                         level--;
 766                     }
 767                     i++;
 768                 }
 769             }
 770             atoms->nres = atoms->nr;
 771
 772             if (i != atoms->nr)
 773             {
 774                 gmx_fatal(FARGS, "Internal inconsistency in Espresso routines, read %d atoms, expected %d atoms", i, atoms->nr);
 775             }
 776         }
 777         else if (level == 1 && std::strcmp(word, "variable") == 0 && !bFoundVariable)
 778         {
 779             bFoundVariable = TRUE;
 780             level         += check_open_parenthesis(fp, r, infile, "variable");
 781             while (level == 2 && (r = get_espresso_word(fp, word)))
 782             {
 783                 if (level == 2 && std::strcmp(word, "box_l") == 0)
 784                 {
 785                     for (m = 0; m < 3; m++)
 786                     {
 787                         r = get_espresso_word(fp, word);
 788                         sscanf(word, "%lf", &d);
 789                         box[m][m] = d;
 790                     }
 791                     level += check_close_parenthesis(fp, r, infile, "box_l");
 792                 }
 793             }
 794         }
 795         else if (r == 2)
 796         {
 797             level++;
 798         }
 799         else if (r == 3)
 800         {
 801             level--;
 802         }
 803     }
 804
 805     if (!bFoundParticles)
 806     {
 807         fprintf(stderr, "Did not find a particles section in Espresso file '%s'\n",
 808                 infile);
 809     }
 810
 811     gmx_fio_fclose(fp);
 812 }
 813
 814 static int get_espresso_coordnum(const char *infile)
 815 {
 816     FILE    *fp;
 817     char     word[STRLEN];
 818     int      natoms, level, r;
 819     gmx_bool bFoundParticles;
 820
 821     natoms = 0;
 822
 823     fp = gmx_fio_fopen(infile, "r");
 824
 825     bFoundParticles = FALSE;
 826     level           = 0;
 827     while ((r = get_espresso_word(fp, word)) && !bFoundParticles)
 828     {
 829         if (level == 1 && strcmp(word, "particles") == 0 && !bFoundParticles)
 830         {
 831             bFoundParticles = TRUE;
 832             level          += check_open_parenthesis(fp, r, infile, "particles");
 833             while (level > 0 && (r = get_espresso_word(fp, word)))
 834             {
 835                 if (r == 2)
 836                 {
 837                     level++;
 838                     if (level == 2)
 839                     {
 840                         natoms++;
 841                     }
 842                 }
 843                 else if (r == 3)
 844                 {
 845                     level--;
 846                 }
 847             }
 848         }
 849         else if (r == 2)
 850         {
 851             level++;
 852         }
 853         else if (r == 3)
 854         {
 855             level--;
 856         }
 857     }
 858     if (!bFoundParticles)
 859     {
 860         fprintf(stderr, "Did not find a particles section in Espresso file '%s'\n",
 861                 infile);
 862     }
 863
 864     gmx_fio_fclose(fp);
 865
 866     return natoms;
 867 }
 868
 869 static void write_espresso_conf_indexed(FILE *out, const char *title,
 870                                         t_atoms *atoms, int nx, atom_id *index,
 871                                         rvec *x, rvec *v, matrix box)
 872 {
 873     int i, j;
 874
 875     fprintf(out, "# %s\n", title);
 876     if (TRICLINIC(box))
 877     {
 878         gmx_warning("The Espresso format does not support triclinic unit-cells");
 879     }
 880     fprintf(out, "{variable {box_l %f %f %f}}\n", box[0][0], box[1][1], box[2][2]);
 881
 882     fprintf(out, "{particles {id pos type q%s}\n", v ? " v" : "");
 883     for (i = 0; i < nx; i++)
 884     {
 885         if (index)
 886         {
 887             j = index[i];
 888         }
 889         else
 890         {
 891             j = i;
 892         }
 893         fprintf(out, "\t{%d %f %f %f %d %g",
 894                 j, x[j][XX], x[j][YY], x[j][ZZ],
 895                 atoms->atom[j].type, atoms->atom[j].q);
 896         if (v)
 897         {
 898             fprintf(out, " %f %f %f", v[j][XX], v[j][YY], v[j][ZZ]);
 899         }
 900         fprintf(out, "}\n");
 901     }
 902     fprintf(out, "}\n");
 903 }
 904
 905 static void get_coordnum_fp (FILE *in, char *title, int *natoms)
 906 {
 907     char line[STRLEN+1];
 908
 909     fgets2 (title, STRLEN, in);
 910     fgets2 (line, STRLEN, in);
 911     if (sscanf (line, "%d", natoms) != 1)
 912     {
 913         gmx_fatal(FARGS, "gro file does not have the number of atoms on the second line");
 914     }
 915 }
 916
 917 static void get_coordnum (const char *infile, int *natoms)
 918 {
 919     FILE *in;
 920     char  title[STRLEN];
 921
 922     in = gmx_fio_fopen(infile, "r");
 923     get_coordnum_fp(in, title, natoms);
 924     gmx_fio_fclose (in);
 925 }
 926
 927 static gmx_bool get_w_conf(FILE *in, const char *infile, char *title,
 928                            t_symtab *symtab, t_atoms *atoms, int *ndec,
 929                            rvec x[], rvec *v, matrix box)
 930 {
 931     char       name[6];
 932     char       resname[6], oldresname[6];
 933     char       line[STRLEN+1], *ptr;
 934     char       buf[256];
 935     double     x1, y1, z1, x2, y2, z2;
 936     rvec       xmin, xmax;
 937     int        natoms, i, m, resnr, newres, oldres, ddist, c;
 938     gmx_bool   bFirst, bVel;
 939     char      *p1, *p2, *p3;
 940
 941     newres  = -1;
 942     oldres  = NOTSET; /* Unlikely number for the first residue! */
 943     ddist   = 0;
 944
 945     /* Read the title and number of atoms */
 946     get_coordnum_fp(in, title, &natoms);
 947
 948     if (natoms > atoms->nr)
 949     {
 950         gmx_fatal(FARGS, "gro file contains more atoms (%d) than expected (%d)",
 951                   natoms, atoms->nr);
 952     }
 953     else if (natoms <  atoms->nr)
 954     {
 955         fprintf(stderr, "Warning: gro file contains less atoms (%d) than expected"
 956                 " (%d)\n", natoms, atoms->nr);
 957     }
 958
 959     bFirst = TRUE;
 960
 961     bVel = FALSE;
 962
 963     resname[0]     = '\0';
 964     oldresname[0]  = '\0';
 965
 966     /* just pray the arrays are big enough */
 967     for (i = 0; (i < natoms); i++)
 968     {
 969         if ((fgets2 (line, STRLEN, in)) == NULL)
 970         {
 971             gmx_fatal(FARGS, "Unexpected end of file in file %s at line %d",
 972                       infile, i+2);
 973         }
 974         if (strlen(line) < 39)
 975         {
 976             gmx_fatal(FARGS, "Invalid line in %s for atom %d:\n%s", infile, i+1, line);
 977         }
 978
 979         /* determine read precision from distance between periods
 980            (decimal points) */
 981         if (bFirst)
 982         {
 983             bFirst = FALSE;
 984             p1     = strchr(line, '.');
 985             if (p1 == NULL)
 986             {
 987                 gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
 988             }
 989             p2 = strchr(&p1[1], '.');
 990             if (p2 == NULL)
 991             {
 992                 gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
 993             }
 994             ddist = p2 - p1;
 995             *ndec = ddist - 5;
 996
 997             p3 = strchr(&p2[1], '.');
 998             if (p3 == NULL)
 999             {
1000                 gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
1001             }
1002
1003             if (p3 - p2 != ddist)
1004             {
1005                 gmx_fatal(FARGS, "The spacing of the decimal points in file %s is not consistent for x, y and z", infile);
1006             }
1007         }
1008
1009         /* residue number*/
1010         memcpy(name, line, 5);
1011         name[5] = '\0';
1012         sscanf(name, "%d", &resnr);
1013         sscanf(line+5, "%5s", resname);
1014
1015         if (resnr != oldres || strncmp(resname, oldresname, sizeof(resname)))
1016         {
1017             oldres = resnr;
1018             newres++;
1019             if (newres >= natoms)
1020             {
1021                 gmx_fatal(FARGS, "More residues than atoms in %s (natoms = %d)",
1022                           infile, natoms);
1023             }
1024             atoms->atom[i].resind = newres;
1025             t_atoms_set_resinfo(atoms, i, symtab, resname, resnr, ' ', 0, ' ');
1026         }
1027         else
1028         {
1029             atoms->atom[i].resind = newres;
1030         }
1031
1032         /* atomname */
1033         std::memcpy(name, line+10, 5);
1034         atoms->atomname[i] = put_symtab(symtab, name);
1035
1036         /* Copy resname to oldresname after we are done with the sanity check above */
1037         std::strncpy(oldresname, resname, sizeof(oldresname));
1038
1039         /* eventueel controle atomnumber met i+1 */
1040
1041         /* coordinates (start after residue data) */
1042         ptr = line + 20;
1043         /* Read fixed format */
1044         for (m = 0; m < DIM; m++)
1045         {
1046             for (c = 0; (c < ddist && ptr[0]); c++)
1047             {
1048                 buf[c] = ptr[0];
1049                 ptr++;
1050             }
1051             buf[c] = '\0';
1052             if (sscanf (buf, "%lf %lf", &x1, &x2) != 1)
1053             {
1054                 gmx_fatal(FARGS, "Something is wrong in the coordinate formatting of file %s. Note that gro is fixed format (see the manual)", infile);
1055             }
1056             else
1057             {
1058                 x[i][m] = x1;
1059             }
1060         }
1061
1062         /* velocities (start after residues and coordinates) */
1063         if (v)
1064         {
1065             /* Read fixed format */
1066             for (m = 0; m < DIM; m++)
1067             {
1068                 for (c = 0; (c < ddist && ptr[0]); c++)
1069                 {
1070                     buf[c] = ptr[0];
1071                     ptr++;
1072                 }
1073                 buf[c] = '\0';
1074                 if (sscanf (buf, "%lf", &x1) != 1)
1075                 {
1076                     v[i][m] = 0;
1077                 }
1078                 else
1079                 {
1080                     v[i][m] = x1;
1081                     bVel    = TRUE;
1082                 }
1083             }
1084         }
1085     }
1086     atoms->nres = newres + 1;
1087
1088     /* box */
1089     fgets2 (line, STRLEN, in);
1090     if (sscanf (line, "%lf%lf%lf", &x1, &y1, &z1) != 3)
1091     {
1092         gmx_warning("Bad box in file %s", infile);
1093
1094         /* Generate a cubic box */
1095         for (m = 0; (m < DIM); m++)
1096         {
1097             xmin[m] = xmax[m] = x[0][m];
1098         }
1099         for (i = 1; (i < atoms->nr); i++)
1100         {
1101             for (m = 0; (m < DIM); m++)
1102             {
1103                 xmin[m] = std::min(xmin[m], x[i][m]);
1104                 xmax[m] = std::max(xmax[m], x[i][m]);
1105             }
1106         }
1107         for (i = 0; i < DIM; i++)
1108         {
1109             for (m = 0; m < DIM; m++)
1110             {
1111                 box[i][m] = 0.0;
1112             }
1113         }
1114         for (m = 0; (m < DIM); m++)
1115         {
1116             box[m][m] = (xmax[m]-xmin[m]);
1117         }
1118         fprintf(stderr, "Generated a cubic box %8.3f x %8.3f x %8.3f\n",
1119                 box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
1120     }
1121     else
1122     {
1123         /* We found the first three values, the diagonal elements */
1124         box[XX][XX] = x1;
1125         box[YY][YY] = y1;
1126         box[ZZ][ZZ] = z1;
1127         if (sscanf (line, "%*f%*f%*f%lf%lf%lf%lf%lf%lf",
1128                     &x1, &y1, &z1, &x2, &y2, &z2) != 6)
1129         {
1130             x1 = y1 = z1 = x2 = y2 = z2 = 0.0;
1131         }
1132         box[XX][YY] = x1;
1133         box[XX][ZZ] = y1;
1134         box[YY][XX] = z1;
1135         box[YY][ZZ] = x2;
1136         box[ZZ][XX] = y2;
1137         box[ZZ][YY] = z2;
1138     }
1139
1140     return bVel;
1141 }
1142
1143 static void read_whole_conf(const char *infile, char *title,
1144                             t_atoms *atoms, rvec x[], rvec *v, matrix box)
1145 {
1146     FILE    *in;
1147     int      ndec;
1148     t_symtab symtab;
1149
1150     /* open file */
1151     in = gmx_fio_fopen(infile, "r");
1152
1153     open_symtab(&symtab);
1154     get_w_conf(in, infile, title, &symtab, atoms, &ndec, x, v, box);
1155     /* We can't free the symbols, as they are still used in atoms, so
1156      * the only choice is to leak them. */
1157     free_symtab(&symtab);
1158
1159     gmx_fio_fclose(in);
1160 }
1161
1162 static gmx_bool gmx_one_before_eof(FILE *fp)
1163 {
1164     char     data[4];
1165     gmx_bool beof;
1166
1167     if ((beof = fread(data, 1, 1, fp)) == 1)
1168     {
1169         gmx_fseek(fp, -1, SEEK_CUR);
1170     }
1171     return !beof;
1172 }
1173
1174 gmx_bool gro_next_x_or_v(FILE *status, t_trxframe *fr)
1175 {
1176     t_atoms  atoms;
1177     t_symtab symtab;
1178     char     title[STRLEN], *p;
1179     double   tt;
1180     int      ndec = 0, i;
1181
1182     if (gmx_one_before_eof(status))
1183     {
1184         return FALSE;
1185     }
1186
1187     open_symtab(&symtab);
1188     atoms.nr = fr->natoms;
1189     snew(atoms.atom, fr->natoms);
1190     atoms.nres = fr->natoms;
1191     snew(atoms.resinfo, fr->natoms);
1192     snew(atoms.atomname, fr->natoms);
1193
1194     fr->bV    = get_w_conf(status, title, title, &symtab, &atoms, &ndec, fr->x, fr->v, fr->box);
1195     fr->bPrec = TRUE;
1196     fr->prec  = 1;
1197     /* prec = 10^ndec: */
1198     for (i = 0; i < ndec; i++)
1199     {
1200         fr->prec *= 10;
1201     }
1202     fr->title  = title;
1203     fr->bTitle = TRUE;
1204     fr->bX     = TRUE;
1205     fr->bBox   = TRUE;
1206
1207     sfree(atoms.atom);
1208     sfree(atoms.resinfo);
1209     sfree(atoms.atomname);
1210     done_symtab(&symtab);
1211
1212     if ((p = strstr(title, "t=")) != NULL)
1213     {
1214         p += 2;
1215         if (sscanf(p, "%lf", &tt) == 1)
1216         {
1217             fr->time  = tt;
1218             fr->bTime = TRUE;
1219         }
1220         else
1221         {
1222             fr->time  = 0;
1223             fr->bTime = FALSE;
1224         }
1225     }
1226
1227     if (atoms.nr != fr->natoms)
1228     {
1229         gmx_fatal(FARGS, "Number of atoms in gro frame (%d) doesn't match the number in the previous frame (%d)", atoms.nr, fr->natoms);
1230     }
1231
1232     return TRUE;
1233 }
1234
1235 int gro_first_x_or_v(FILE *status, t_trxframe *fr)
1236 {
1237     char title[STRLEN];
1238
1239     frewind(status);
1240     fprintf(stderr, "Reading frames from gro file");
1241     get_coordnum_fp(status, title, &fr->natoms);
1242     frewind(status);
1243     fprintf(stderr, " '%s', %d atoms.\n", title, fr->natoms);
1244     fr->bTitle = TRUE;
1245     fr->title  = title;
1246     if (fr->natoms == 0)
1247     {
1248         gmx_file("No coordinates in gro file");
1249     }
1250
1251     snew(fr->x, fr->natoms);
1252     snew(fr->v, fr->natoms);
1253     gro_next_x_or_v(status, fr);
1254
1255     return fr->natoms;
1256 }
1257
1258 static void make_hconf_format(int pr, gmx_bool bVel, char format[])
1259 {
1260     int l, vpr;
1261
1262     /* build format string for printing,
1263        something like "%8.3f" for x and "%8.4f" for v */
1264     if (pr < 0)
1265     {
1266         pr = 0;
1267     }
1268     if (pr > 30)
1269     {
1270         pr = 30;
1271     }
1272     l   = pr+5;
1273     vpr = pr+1;
1274     if (bVel)
1275     {
1276         sprintf(format, "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
1277                 l, pr, l, pr, l, pr, l, vpr, l, vpr, l, vpr);
1278     }
1279     else
1280     {
1281         sprintf(format, "%%%d.%df%%%d.%df%%%d.%df\n", l, pr, l, pr, l, pr);
1282     }
1283
1284 }
1285
1286 static void write_hconf_box(FILE *out, int pr, matrix box)
1287 {
1288     char format[100];
1289     int  l;
1290
1291     if (pr < 5)
1292     {
1293         pr = 5;
1294     }
1295     l = pr+5;
1296
1297     if (box[XX][YY] || box[XX][ZZ] || box[YY][XX] || box[YY][ZZ] ||
1298         box[ZZ][XX] || box[ZZ][YY])
1299     {
1300         sprintf(format, "%%%d.%df%%%d.%df%%%d.%df"
1301                 "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
1302                 l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr);
1303         fprintf(out, format,
1304                 box[XX][XX], box[YY][YY], box[ZZ][ZZ],
1305                 box[XX][YY], box[XX][ZZ], box[YY][XX],
1306                 box[YY][ZZ], box[ZZ][XX], box[ZZ][YY]);
1307     }
1308     else
1309     {
1310         sprintf(format, "%%%d.%df%%%d.%df%%%d.%df\n", l, pr, l, pr, l, pr);
1311         fprintf(out, format,
1312                 box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
1313     }
1314 }
1315
1316 void write_hconf_indexed_p(FILE *out, const char *title, t_atoms *atoms,
1317                            int nx, const atom_id index[], int pr,
1318                            rvec *x, rvec *v, matrix box)
1319 {
1320     char resnm[6], nm[6], format[100];
1321     int  ai, i, resind, resnr;
1322
1323     bromacs(format, 99);
1324     fprintf (out, "%s\n", (title && title[0]) ? title : format);
1325     fprintf (out, "%5d\n", nx);
1326
1327     make_hconf_format(pr, v != NULL, format);
1328
1329     for (i = 0; (i < nx); i++)
1330     {
1331         ai = index[i];
1332
1333         resind = atoms->atom[ai].resind;
1334         std::strncpy(resnm, " ??? ", sizeof(resnm)-1);
1335         if (resind < atoms->nres)
1336         {
1337             std::strncpy(resnm, *atoms->resinfo[resind].name, sizeof(resnm)-1);
1338             resnr = atoms->resinfo[resind].nr;
1339         }
1340         else
1341         {
1342             std::strncpy(resnm, " ??? ", sizeof(resnm)-1);
1343             resnr = resind + 1;
1344         }
1345
1346         if (atoms->atom)
1347         {
1348             std::strncpy(nm, *atoms->atomname[ai], sizeof(nm)-1);
1349         }
1350         else
1351         {
1352             std::strncpy(nm, " ??? ", sizeof(nm)-1);
1353         }
1354
1355         fprintf(out, "%5d%-5.5s%5.5s%5d", resnr%100000, resnm, nm, (ai+1)%100000);
1356         /* next fprintf uses built format string */
1357         if (v)
1358         {
1359             fprintf(out, format,
1360                     x[ai][XX], x[ai][YY], x[ai][ZZ], v[ai][XX], v[ai][YY], v[ai][ZZ]);
1361         }
1362         else
1363         {
1364             fprintf(out, format,
1365                     x[ai][XX], x[ai][YY], x[ai][ZZ]);
1366         }
1367     }
1368
1369     write_hconf_box(out, pr, box);
1370
1371     fflush(out);
1372 }
1373
1374 static void write_hconf_mtop(FILE *out, const char *title, gmx_mtop_t *mtop,
1375                              int pr,
1376                              rvec *x, rvec *v, matrix box)
1377 {
1378     char                    format[100];
1379     int                     i, resnr;
1380     gmx_mtop_atomloop_all_t aloop;
1381     t_atom                 *atom;
1382     char                   *atomname, *resname;
1383
1384     bromacs(format, 99);
1385     fprintf (out, "%s\n", (title && title[0]) ? title : format);
1386     fprintf (out, "%5d\n", mtop->natoms);
1387
1388     make_hconf_format(pr, v != NULL, format);
1389
1390     aloop = gmx_mtop_atomloop_all_init(mtop);
1391     while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
1392     {
1393         gmx_mtop_atomloop_all_names(aloop, &atomname, &resnr, &resname);
1394
1395         fprintf(out, "%5d%-5.5s%5.5s%5d",
1396                 resnr%100000, resname, atomname, (i+1)%100000);
1397         /* next fprintf uses built format string */
1398         if (v)
1399         {
1400             fprintf(out, format,
1401                     x[i][XX], x[i][YY], x[i][ZZ], v[i][XX], v[i][YY], v[i][ZZ]);
1402         }
1403         else
1404         {
1405             fprintf(out, format,
1406                     x[i][XX], x[i][YY], x[i][ZZ]);
1407         }
1408     }
1409
1410     write_hconf_box(out, pr, box);
1411
1412     fflush(out);
1413 }
1414
1415 void write_hconf_p(FILE *out, const char *title, t_atoms *atoms, int pr,
1416                    rvec *x, rvec *v, matrix box)
1417 {
1418     atom_id *aa;
1419     int      i;
1420
1421     snew(aa, atoms->nr);
1422     for (i = 0; (i < atoms->nr); i++)
1423     {
1424         aa[i] = i;
1425     }
1426     write_hconf_indexed_p(out, title, atoms, atoms->nr, aa, pr, x, v, box);
1427     sfree(aa);
1428 }
1429
1430 void write_conf_p(const char *outfile, const char *title,
1431                   t_atoms *atoms, int pr,
1432                   rvec *x, rvec *v, matrix box)
1433 {
1434     FILE *out;
1435
1436     out = gmx_fio_fopen(outfile, "w");
1437     write_hconf_p(out, title, atoms, pr, x, v, box);
1438
1439     gmx_fio_fclose (out);
1440 }
1441
1442 static void write_conf(const char *outfile, const char *title, t_atoms *atoms,
1443                        rvec *x, rvec *v, matrix box)
1444 {
1445     write_conf_p(outfile, title, atoms, 3, x, v, box);
1446 }
1447
1448 void write_sto_conf_indexed(const char *outfile, const char *title,
1449                             t_atoms *atoms,
1450                             rvec x[], rvec *v, int ePBC, matrix box,
1451                             atom_id nindex, atom_id index[])
1452 {
1453     FILE       *out;
1454     int         ftp;
1455     t_trxframe  fr;
1456
1457     ftp = fn2ftp(outfile);
1458     switch (ftp)
1459     {
1460         case efGRO:
1461             out = gmx_fio_fopen(outfile, "w");
1462             write_hconf_indexed_p(out, title, atoms, nindex, index, 3, x, v, box);
1463             gmx_fio_fclose(out);
1464             break;
1465         case efG96:
1466             clear_trxframe(&fr, TRUE);
1467             fr.bTitle = TRUE;
1468             fr.title  = title;
1469             fr.natoms = atoms->nr;
1470             fr.bAtoms = TRUE;
1471             fr.atoms  = atoms;
1472             fr.bX     = TRUE;
1473             fr.x      = x;
1474             if (v)
1475             {
1476                 fr.bV = TRUE;
1477                 fr.v  = v;
1478             }
1479             fr.bBox = TRUE;
1480             copy_mat(box, fr.box);
1481             out = gmx_fio_fopen(outfile, "w");
1482             write_g96_conf(out, &fr, nindex, index);
1483             gmx_fio_fclose(out);
1484             break;
1485         case efPDB:
1486         case efBRK:
1487         case efENT:
1488         case efPQR:
1489             out = gmx_fio_fopen(outfile, "w");
1490             write_pdbfile_indexed(out, title, atoms, x, ePBC, box, ' ', -1, nindex, index, NULL, TRUE);
1491             gmx_fio_fclose(out);
1492             break;
1493         case efESP:
1494             out = gmx_fio_fopen(outfile, "w");
1495             write_espresso_conf_indexed(out, title, atoms, nindex, index, x, v, box);
1496             gmx_fio_fclose(out);
1497             break;
1498         case efTPR:
1499             gmx_fatal(FARGS, "Sorry, can not write a topology to %s", outfile);
1500             break;
1501         default:
1502             gmx_incons("Not supported in write_sto_conf_indexed");
1503     }
1504 }
1505
1506 void write_sto_conf(const char *outfile, const char *title, t_atoms *atoms,
1507                     rvec x[], rvec *v, int ePBC, matrix box)
1508 {
1509     FILE       *out;
1510     int         ftp;
1511     t_trxframe  fr;
1512
1513     ftp = fn2ftp(outfile);
1514     switch (ftp)
1515     {
1516         case efGRO:
1517             write_conf(outfile, title, atoms, x, v, box);
1518             break;
1519         case efG96:
1520             clear_trxframe(&fr, TRUE);
1521             fr.bTitle = TRUE;
1522             fr.title  = title;
1523             fr.natoms = atoms->nr;
1524             fr.bAtoms = TRUE;
1525             fr.atoms  = atoms;
1526             fr.bX     = TRUE;
1527             fr.x      = x;
1528             if (v)
1529             {
1530                 fr.bV = TRUE;
1531                 fr.v  = v;
1532             }
1533             fr.bBox = TRUE;
1534             copy_mat(box, fr.box);
1535             out = gmx_fio_fopen(outfile, "w");
1536             write_g96_conf(out, &fr, -1, NULL);
1537             gmx_fio_fclose(out);
1538             break;
1539         case efPDB:
1540         case efBRK:
1541         case efENT:
1542             out = gmx_fio_fopen(outfile, "w");
1543             write_pdbfile(out, title, atoms, x, ePBC, box, ' ', -1, NULL, TRUE);
1544             gmx_fio_fclose(out);
1545             break;
1546         case efESP:
1547             out = gmx_fio_fopen(outfile, "w");
1548             write_espresso_conf_indexed(out, title, atoms, atoms->nr, NULL, x, v, box);
1549             gmx_fio_fclose(out);
1550             break;
1551         case efTPR:
1552             gmx_fatal(FARGS, "Sorry, can not write a topology to %s", outfile);
1553             break;
1554         default:
1555             gmx_incons("Not supported in write_sto_conf");
1556     }
1557 }
1558
1559 void write_sto_conf_mtop(const char *outfile, const char *title,
1560                          gmx_mtop_t *mtop,
1561                          rvec x[], rvec *v, int ePBC, matrix box)
1562 {
1563     int     ftp;
1564     FILE   *out;
1565     t_atoms atoms;
1566
1567     ftp = fn2ftp(outfile);
1568     switch (ftp)
1569     {
1570         case efGRO:
1571             out = gmx_fio_fopen(outfile, "w");
1572             write_hconf_mtop(out, title, mtop, 3, x, v, box);
1573             gmx_fio_fclose(out);
1574             break;
1575         default:
1576             /* This is a brute force approach which requires a lot of memory.
1577              * We should implement mtop versions of all writing routines.
1578              */
1579             atoms = gmx_mtop_global_atoms(mtop);
1580
1581             write_sto_conf(outfile, title, &atoms, x, v, ePBC, box);
1582
1583             done_atom(&atoms);
1584             break;
1585     }
1586 }
1587
1588 void get_stx_coordnum(const char *infile, int *natoms)
1589 {
1590     FILE      *in;
1591     int        ftp, tpxver, tpxgen;
1592     t_trxframe fr;
1593     char       g96_line[STRLEN+1];
1594
1595     ftp = fn2ftp(infile);
1596     range_check(ftp, 0, efNR);
1597     switch (ftp)
1598     {
1599         case efGRO:
1600             get_coordnum(infile, natoms);
1601             break;
1602         case efG96:
1603             in        = gmx_fio_fopen(infile, "r");
1604             fr.title  = NULL;
1605             fr.natoms = -1;
1606             fr.atoms  = NULL;
1607             fr.x      = NULL;
1608             fr.v      = NULL;
1609             fr.f      = NULL;
1610             *natoms   = read_g96_conf(in, infile, &fr, g96_line);
1611             gmx_fio_fclose(in);
1612             break;
1613         case efPDB:
1614         case efBRK:
1615         case efENT:
1616             in = gmx_fio_fopen(infile, "r");
1617             get_pdb_coordnum(in, natoms);
1618             gmx_fio_fclose(in);
1619             break;
1620         case efESP:
1621             *natoms = get_espresso_coordnum(infile);
1622             break;
1623         case efTPR:
1624         {
1625             t_tpxheader tpx;
1626
1627             read_tpxheader(infile, &tpx, TRUE, &tpxver, &tpxgen);
1628             *natoms = tpx.natoms;
1629             break;
1630         }
1631         default:
1632             gmx_fatal(FARGS, "File type %s not supported in get_stx_coordnum",
1633                       ftp2ext(ftp));
1634     }
1635 }
1636
1637 void read_stx_conf(const char *infile, char *title, t_atoms *atoms,
1638                    rvec x[], rvec *v, int *ePBC, matrix box)
1639 {
1640     FILE       *in;
1641     gmx_mtop_t *mtop;
1642     t_topology  top;
1643     t_trxframe  fr;
1644     int         i, ftp, natoms;
1645     char        g96_line[STRLEN+1];
1646
1647     if (atoms->nr == 0)
1648     {
1649         fprintf(stderr, "Warning: Number of atoms in %s is 0\n", infile);
1650     }
1651     else if (atoms->atom == NULL)
1652     {
1653         gmx_mem("Uninitialized array atom");
1654     }
1655
1656     if (ePBC)
1657     {
1658         *ePBC = -1;
1659     }
1660
1661     ftp = fn2ftp(infile);
1662     switch (ftp)
1663     {
1664         case efGRO:
1665             read_whole_conf(infile, title, atoms, x, v, box);
1666             break;
1667         case efG96:
1668             fr.title  = NULL;
1669             fr.natoms = atoms->nr;
1670             fr.atoms  = atoms;
1671             fr.x      = x;
1672             fr.v      = v;
1673             fr.f      = NULL;
1674             in        = gmx_fio_fopen(infile, "r");
1675             read_g96_conf(in, infile, &fr, g96_line);
1676             gmx_fio_fclose(in);
1677             copy_mat(fr.box, box);
1678             std::strncpy(title, fr.title, STRLEN);
1679             break;
1680         case efPDB:
1681         case efBRK:
1682         case efENT:
1683             read_pdb_conf(infile, title, atoms, x, ePBC, box, TRUE, NULL);
1684             break;
1685         case efESP:
1686             read_espresso_conf(infile, title, atoms, x, v, box);
1687             break;
1688         case efTPR:
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 }