src/gmxlib/confio.c

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