src/gromacs/fileio/confio.c

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