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