src/gromacs/fileio/confio.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "confio.h"
  40
  41 #include <cerrno>
  42 #include <cmath>
  43 #include <cstdio>
  44 #include <cstdlib>
  45 #include <cstring>
  46
  47 #include <algorithm>
  48
  49 #include "gromacs/fileio/filenm.h"
  50 #include "gromacs/fileio/g96io.h"
  51 #include "gromacs/fileio/gmxfio.h"
  52 #include "gromacs/fileio/pdbio.h"
  53 #include "gromacs/fileio/tpxio.h"
  54 #include "gromacs/fileio/trx.h"
  55 #include "gromacs/fileio/trxio.h"
  56 #include "gromacs/legacyheaders/copyrite.h"
  57 #include "gromacs/math/vec.h"
  58 #include "gromacs/pbcutil/pbc.h"
  59 #include "gromacs/topology/atomprop.h"
  60 #include "gromacs/topology/atoms.h"
  61 #include "gromacs/topology/block.h"
  62 #include "gromacs/topology/mtop_util.h"
  63 #include "gromacs/topology/symtab.h"
  64 #include "gromacs/topology/topology.h"
  65 #include "gromacs/utility/cstringutil.h"
  66 #include "gromacs/utility/fatalerror.h"
  67 #include "gromacs/utility/smalloc.h"
  68
  69 static int get_espresso_word(FILE *fp, char word[])
  70 {
  71     int  ret, nc, i;
  72
  73     ret = 0;
  74     nc  = 0;
  75
  76     do
  77     {
  78         i = fgetc(fp);
  79         if (i != EOF)
  80         {
  81             if (i == ' ' || i == '\n' || i == '\t')
  82             {
  83                 if (nc > 0)
  84                 {
  85                     ret = 1;
  86                 }
  87             }
  88             else if (i == '{')
  89             {
  90                 if (nc == 0)
  91                 {
  92                     word[nc++] = '{';
  93                 }
  94                 ret = 2;
  95             }
  96             else if (i == '}')
  97             {
  98                 if (nc == 0)
  99                 {
 100                     word[nc++] = '}';
 101                 }
 102                 ret = 3;
 103             }
 104             else
 105             {
 106                 word[nc++] = (char)i;
 107             }
 108         }
 109     }
 110     while (i != EOF && ret == 0);
 111
 112     word[nc] = '\0';
 113
 114     /*  printf("'%s'\n",word); */
 115
 116     return ret;
 117 }
 118
 119 static int check_open_parenthesis(FILE *fp, int r,
 120                                   const char *infile, const char *keyword)
 121 {
 122     int  level_inc;
 123     char word[STRLEN];
 124
 125     level_inc = 0;
 126     if (r == 2)
 127     {
 128         level_inc++;
 129     }
 130     else
 131     {
 132         r = get_espresso_word(fp, word);
 133         if (r == 2)
 134         {
 135             level_inc++;
 136         }
 137         else
 138         {
 139             gmx_fatal(FARGS, "Expected '{' after '%s' in file '%s'",
 140                       keyword, infile);
 141         }
 142     }
 143
 144     return level_inc;
 145 }
 146
 147 static int check_close_parenthesis(FILE *fp, int r,
 148                                    const char *infile, const char *keyword)
 149 {
 150     int  level_inc;
 151     char word[STRLEN];
 152
 153     level_inc = 0;
 154     if (r == 3)
 155     {
 156         level_inc--;
 157     }
 158     else
 159     {
 160         r = get_espresso_word(fp, word);
 161         if (r == 3)
 162         {
 163             level_inc--;
 164         }
 165         else
 166         {
 167             gmx_fatal(FARGS, "Expected '}' after section '%s' in file '%s'",
 168                       keyword, infile);
 169         }
 170     }
 171
 172     return level_inc;
 173 }
 174
 175 enum {
 176     espID, espPOS, espTYPE, espQ, espV, espF, espMOLECULE, espNR
 177 };
 178 const char *esp_prop[espNR] = {
 179     "id", "pos", "type", "q", "v", "f",
 180     "molecule"
 181 };
 182
 183 static void read_espresso_conf(const char *infile, char *title,
 184                                t_atoms *atoms, rvec x[], rvec *v, matrix box)
 185 {
 186     t_symtab *symtab = NULL;
 187     FILE     *fp;
 188     char      word[STRLEN], buf[STRLEN];
 189     int       level, r, nprop, p, i, m, molnr;
 190     int       prop[32];
 191     double    d;
 192     gmx_bool  bFoundParticles, bFoundProp, bFoundVariable, bMol;
 193
 194     if (!symtab)
 195     {
 196         snew(symtab, 1);
 197         open_symtab(symtab);
 198     }
 199     // TODO: The code does not understand titles it writes...
 200     title[0] = '\0';
 201
 202     clear_mat(box);
 203
 204     fp = gmx_fio_fopen(infile, "r");
 205
 206     bFoundParticles = FALSE;
 207     bFoundVariable  = FALSE;
 208     bMol            = FALSE;
 209     level           = 0;
 210     while ((r = get_espresso_word(fp, word)))
 211     {
 212         if (level == 1 && std::strcmp(word, "particles") == 0 && !bFoundParticles)
 213         {
 214             bFoundParticles = TRUE;
 215             level          += check_open_parenthesis(fp, r, infile, "particles");
 216             nprop           = 0;
 217             while (level == 2 && (r = get_espresso_word(fp, word)))
 218             {
 219                 bFoundProp = FALSE;
 220                 for (p = 0; p < espNR; p++)
 221                 {
 222                     if (strcmp(word, esp_prop[p]) == 0)
 223                     {
 224                         bFoundProp    = TRUE;
 225                         prop[nprop++] = p;
 226                         /* printf("  prop[%d] = %s\n",nprop-1,esp_prop[prop[nprop-1]]); */
 227                     }
 228                 }
 229                 if (!bFoundProp && word[0] != '}')
 230                 {
 231                     gmx_fatal(FARGS, "Can not read Espresso files with particle property '%s'", word);
 232                 }
 233                 if (bFoundProp && p == espMOLECULE)
 234                 {
 235                     bMol = TRUE;
 236                 }
 237                 if (r == 3)
 238                 {
 239                     level--;
 240                 }
 241             }
 242
 243             i = 0;
 244             while (level > 0 && (r = get_espresso_word(fp, word)))
 245             {
 246                 if (r == 2)
 247                 {
 248                     level++;
 249                 }
 250                 else if (r == 3)
 251                 {
 252                     level--;
 253                 }
 254                 if (level == 2)
 255                 {
 256                     for (p = 0; p < nprop; p++)
 257                     {
 258                         switch (prop[p])
 259                         {
 260                             case espID:
 261                                 r = get_espresso_word(fp, word);
 262                                 /* Not used */
 263                                 break;
 264                             case espPOS:
 265                                 for (m = 0; m < 3; m++)
 266                                 {
 267                                     r = get_espresso_word(fp, word);
 268                                     sscanf(word, "%lf", &d);
 269                                     x[i][m] = d;
 270                                 }
 271                                 break;
 272                             case espTYPE:
 273                                 r                   = get_espresso_word(fp, word);
 274                                 atoms->atom[i].type = std::strtol(word, NULL, 10);
 275                                 break;
 276                             case espQ:
 277                                 r = get_espresso_word(fp, word);
 278                                 sscanf(word, "%lf", &d);
 279                                 atoms->atom[i].q = d;
 280                                 break;
 281                             case espV:
 282                                 for (m = 0; m < 3; m++)
 283                                 {
 284                                     r = get_espresso_word(fp, word);
 285                                     sscanf(word, "%lf", &d);
 286                                     v[i][m] = d;
 287                                 }
 288                                 break;
 289                             case espF:
 290                                 for (m = 0; m < 3; m++)
 291                                 {
 292                                     r = get_espresso_word(fp, word);
 293                                     /* not used */
 294                                 }
 295                                 break;
 296                             case espMOLECULE:
 297                                 r     = get_espresso_word(fp, word);
 298                                 molnr = std::strtol(word, NULL, 10);
 299                                 if (i == 0 ||
 300                                     atoms->resinfo[atoms->atom[i-1].resind].nr != molnr)
 301                                 {
 302                                     atoms->atom[i].resind =
 303                                         (i == 0 ? 0 : atoms->atom[i-1].resind+1);
 304                                     atoms->resinfo[atoms->atom[i].resind].nr       = molnr;
 305                                     atoms->resinfo[atoms->atom[i].resind].ic       = ' ';
 306                                     atoms->resinfo[atoms->atom[i].resind].chainid  = ' ';
 307                                     atoms->resinfo[atoms->atom[i].resind].chainnum = molnr; /* Not sure if this is right? */
 308                                 }
 309                                 else
 310                                 {
 311                                     atoms->atom[i].resind = atoms->atom[i-1].resind;
 312                                 }
 313                                 break;
 314                         }
 315                     }
 316                     /* Generate an atom name from the particle type */
 317                     sprintf(buf, "T%d", atoms->atom[i].type);
 318                     atoms->atomname[i] = put_symtab(symtab, buf);
 319                     if (bMol)
 320                     {
 321                         if (i == 0 || atoms->atom[i].resind != atoms->atom[i-1].resind)
 322                         {
 323                             atoms->resinfo[atoms->atom[i].resind].name =
 324                                 put_symtab(symtab, "MOL");
 325                         }
 326                     }
 327                     else
 328                     {
 329                         /* Residue number is the atom number */
 330                         atoms->atom[i].resind = i;
 331                         /* Generate an residue name from the particle type */
 332                         if (atoms->atom[i].type < 26)
 333                         {
 334                             sprintf(buf, "T%c", 'A'+atoms->atom[i].type);
 335                         }
 336                         else
 337                         {
 338                             sprintf(buf, "T%c%c",
 339                                     'A'+atoms->atom[i].type/26, 'A'+atoms->atom[i].type%26);
 340                         }
 341                         t_atoms_set_resinfo(atoms, i, symtab, buf, i, ' ', 0, ' ');
 342                     }
 343
 344                     if (r == 3)
 345                     {
 346                         level--;
 347                     }
 348                     i++;
 349                 }
 350             }
 351             atoms->nres = atoms->nr;
 352
 353             if (i != atoms->nr)
 354             {
 355                 gmx_fatal(FARGS, "Internal inconsistency in Espresso routines, read %d atoms, expected %d atoms", i, atoms->nr);
 356             }
 357         }
 358         else if (level == 1 && std::strcmp(word, "variable") == 0 && !bFoundVariable)
 359         {
 360             bFoundVariable = TRUE;
 361             level         += check_open_parenthesis(fp, r, infile, "variable");
 362             while (level == 2 && (r = get_espresso_word(fp, word)))
 363             {
 364                 if (level == 2 && std::strcmp(word, "box_l") == 0)
 365                 {
 366                     for (m = 0; m < 3; m++)
 367                     {
 368                         r = get_espresso_word(fp, word);
 369                         sscanf(word, "%lf", &d);
 370                         box[m][m] = d;
 371                     }
 372                     level += check_close_parenthesis(fp, r, infile, "box_l");
 373                 }
 374             }
 375         }
 376         else if (r == 2)
 377         {
 378             level++;
 379         }
 380         else if (r == 3)
 381         {
 382             level--;
 383         }
 384     }
 385
 386     if (!bFoundParticles)
 387     {
 388         fprintf(stderr, "Did not find a particles section in Espresso file '%s'\n",
 389                 infile);
 390     }
 391
 392     gmx_fio_fclose(fp);
 393 }
 394
 395 static int get_espresso_coordnum(const char *infile)
 396 {
 397     FILE    *fp;
 398     char     word[STRLEN];
 399     int      natoms, level, r;
 400     gmx_bool bFoundParticles;
 401
 402     natoms = 0;
 403
 404     fp = gmx_fio_fopen(infile, "r");
 405
 406     bFoundParticles = FALSE;
 407     level           = 0;
 408     while ((r = get_espresso_word(fp, word)) && !bFoundParticles)
 409     {
 410         if (level == 1 && strcmp(word, "particles") == 0 && !bFoundParticles)
 411         {
 412             bFoundParticles = TRUE;
 413             level          += check_open_parenthesis(fp, r, infile, "particles");
 414             while (level > 0 && (r = get_espresso_word(fp, word)))
 415             {
 416                 if (r == 2)
 417                 {
 418                     level++;
 419                     if (level == 2)
 420                     {
 421                         natoms++;
 422                     }
 423                 }
 424                 else if (r == 3)
 425                 {
 426                     level--;
 427                 }
 428             }
 429         }
 430         else if (r == 2)
 431         {
 432             level++;
 433         }
 434         else if (r == 3)
 435         {
 436             level--;
 437         }
 438     }
 439     if (!bFoundParticles)
 440     {
 441         fprintf(stderr, "Did not find a particles section in Espresso file '%s'\n",
 442                 infile);
 443     }
 444
 445     gmx_fio_fclose(fp);
 446
 447     return natoms;
 448 }
 449
 450 static void write_espresso_conf_indexed(FILE *out, const char *title,
 451                                         t_atoms *atoms, int nx, atom_id *index,
 452                                         rvec *x, rvec *v, matrix box)
 453 {
 454     int i, j;
 455
 456     fprintf(out, "# %s\n", title);
 457     if (TRICLINIC(box))
 458     {
 459         gmx_warning("The Espresso format does not support triclinic unit-cells");
 460     }
 461     fprintf(out, "{variable {box_l %f %f %f}}\n", box[0][0], box[1][1], box[2][2]);
 462
 463     fprintf(out, "{particles {id pos type q%s}\n", v ? " v" : "");
 464     for (i = 0; i < nx; i++)
 465     {
 466         if (index)
 467         {
 468             j = index[i];
 469         }
 470         else
 471         {
 472             j = i;
 473         }
 474         fprintf(out, "\t{%d %f %f %f %d %g",
 475                 j, x[j][XX], x[j][YY], x[j][ZZ],
 476                 atoms->atom[j].type, atoms->atom[j].q);
 477         if (v)
 478         {
 479             fprintf(out, " %f %f %f", v[j][XX], v[j][YY], v[j][ZZ]);
 480         }
 481         fprintf(out, "}\n");
 482     }
 483     fprintf(out, "}\n");
 484 }
 485
 486 static void get_coordnum_fp (FILE *in, char *title, int *natoms)
 487 {
 488     char line[STRLEN+1];
 489
 490     fgets2 (title, STRLEN, in);
 491     fgets2 (line, STRLEN, in);
 492     if (sscanf (line, "%d", natoms) != 1)
 493     {
 494         gmx_fatal(FARGS, "gro file does not have the number of atoms on the second line");
 495     }
 496 }
 497
 498 static void get_coordnum (const char *infile, int *natoms)
 499 {
 500     FILE *in;
 501     char  title[STRLEN];
 502
 503     in = gmx_fio_fopen(infile, "r");
 504     get_coordnum_fp(in, title, natoms);
 505     gmx_fio_fclose (in);
 506 }
 507
 508 static gmx_bool get_w_conf(FILE *in, const char *infile, char *title,
 509                            t_symtab *symtab, t_atoms *atoms, int *ndec,
 510                            rvec x[], rvec *v, matrix box)
 511 {
 512     char       name[6];
 513     char       resname[6], oldresname[6];
 514     char       line[STRLEN+1], *ptr;
 515     char       buf[256];
 516     double     x1, y1, z1, x2, y2, z2;
 517     rvec       xmin, xmax;
 518     int        natoms, i, m, resnr, newres, oldres, ddist, c;
 519     gmx_bool   bFirst, bVel;
 520     char      *p1, *p2, *p3;
 521
 522     newres  = -1;
 523     oldres  = -12345; /* Unlikely number for the first residue! */
 524     ddist   = 0;
 525
 526     /* Read the title and number of atoms */
 527     get_coordnum_fp(in, title, &natoms);
 528
 529     if (natoms > atoms->nr)
 530     {
 531         gmx_fatal(FARGS, "gro file contains more atoms (%d) than expected (%d)",
 532                   natoms, atoms->nr);
 533     }
 534     else if (natoms <  atoms->nr)
 535     {
 536         fprintf(stderr, "Warning: gro file contains less atoms (%d) than expected"
 537                 " (%d)\n", natoms, atoms->nr);
 538     }
 539
 540     bFirst = TRUE;
 541
 542     bVel = FALSE;
 543
 544     resname[0]     = '\0';
 545     oldresname[0]  = '\0';
 546
 547     /* just pray the arrays are big enough */
 548     for (i = 0; (i < natoms); i++)
 549     {
 550         if ((fgets2 (line, STRLEN, in)) == NULL)
 551         {
 552             gmx_fatal(FARGS, "Unexpected end of file in file %s at line %d",
 553                       infile, i+2);
 554         }
 555         if (strlen(line) < 39)
 556         {
 557             gmx_fatal(FARGS, "Invalid line in %s for atom %d:\n%s", infile, i+1, line);
 558         }
 559
 560         /* determine read precision from distance between periods
 561            (decimal points) */
 562         if (bFirst)
 563         {
 564             bFirst = FALSE;
 565             p1     = strchr(line, '.');
 566             if (p1 == NULL)
 567             {
 568                 gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
 569             }
 570             p2 = strchr(&p1[1], '.');
 571             if (p2 == NULL)
 572             {
 573                 gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
 574             }
 575             ddist = p2 - p1;
 576             *ndec = ddist - 5;
 577
 578             p3 = strchr(&p2[1], '.');
 579             if (p3 == NULL)
 580             {
 581                 gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
 582             }
 583
 584             if (p3 - p2 != ddist)
 585             {
 586                 gmx_fatal(FARGS, "The spacing of the decimal points in file %s is not consistent for x, y and z", infile);
 587             }
 588         }
 589
 590         /* residue number*/
 591         memcpy(name, line, 5);
 592         name[5] = '\0';
 593         sscanf(name, "%d", &resnr);
 594         sscanf(line+5, "%5s", resname);
 595
 596         if (resnr != oldres || strncmp(resname, oldresname, sizeof(resname)))
 597         {
 598             oldres = resnr;
 599             newres++;
 600             if (newres >= natoms)
 601             {
 602                 gmx_fatal(FARGS, "More residues than atoms in %s (natoms = %d)",
 603                           infile, natoms);
 604             }
 605             atoms->atom[i].resind = newres;
 606             t_atoms_set_resinfo(atoms, i, symtab, resname, resnr, ' ', 0, ' ');
 607         }
 608         else
 609         {
 610             atoms->atom[i].resind = newres;
 611         }
 612
 613         /* atomname */
 614         std::memcpy(name, line+10, 5);
 615         atoms->atomname[i] = put_symtab(symtab, name);
 616
 617         /* Copy resname to oldresname after we are done with the sanity check above */
 618         std::strncpy(oldresname, resname, sizeof(oldresname));
 619
 620         /* eventueel controle atomnumber met i+1 */
 621
 622         /* coordinates (start after residue data) */
 623         ptr = line + 20;
 624         /* Read fixed format */
 625         for (m = 0; m < DIM; m++)
 626         {
 627             for (c = 0; (c < ddist && ptr[0]); c++)
 628             {
 629                 buf[c] = ptr[0];
 630                 ptr++;
 631             }
 632             buf[c] = '\0';
 633             if (sscanf (buf, "%lf %lf", &x1, &x2) != 1)
 634             {
 635                 gmx_fatal(FARGS, "Something is wrong in the coordinate formatting of file %s. Note that gro is fixed format (see the manual)", infile);
 636             }
 637             else
 638             {
 639                 x[i][m] = x1;
 640             }
 641         }
 642
 643         /* velocities (start after residues and coordinates) */
 644         if (v)
 645         {
 646             /* Read fixed format */
 647             for (m = 0; m < DIM; m++)
 648             {
 649                 for (c = 0; (c < ddist && ptr[0]); c++)
 650                 {
 651                     buf[c] = ptr[0];
 652                     ptr++;
 653                 }
 654                 buf[c] = '\0';
 655                 if (sscanf (buf, "%lf", &x1) != 1)
 656                 {
 657                     v[i][m] = 0;
 658                 }
 659                 else
 660                 {
 661                     v[i][m] = x1;
 662                     bVel    = TRUE;
 663                 }
 664             }
 665         }
 666     }
 667     atoms->nres = newres + 1;
 668
 669     /* box */
 670     fgets2 (line, STRLEN, in);
 671     if (sscanf (line, "%lf%lf%lf", &x1, &y1, &z1) != 3)
 672     {
 673         gmx_warning("Bad box in file %s", infile);
 674
 675         /* Generate a cubic box */
 676         for (m = 0; (m < DIM); m++)
 677         {
 678             xmin[m] = xmax[m] = x[0][m];
 679         }
 680         for (i = 1; (i < atoms->nr); i++)
 681         {
 682             for (m = 0; (m < DIM); m++)
 683             {
 684                 xmin[m] = std::min(xmin[m], x[i][m]);
 685                 xmax[m] = std::max(xmax[m], x[i][m]);
 686             }
 687         }
 688         for (i = 0; i < DIM; i++)
 689         {
 690             for (m = 0; m < DIM; m++)
 691             {
 692                 box[i][m] = 0.0;
 693             }
 694         }
 695         for (m = 0; (m < DIM); m++)
 696         {
 697             box[m][m] = (xmax[m]-xmin[m]);
 698         }
 699         fprintf(stderr, "Generated a cubic box %8.3f x %8.3f x %8.3f\n",
 700                 box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
 701     }
 702     else
 703     {
 704         /* We found the first three values, the diagonal elements */
 705         box[XX][XX] = x1;
 706         box[YY][YY] = y1;
 707         box[ZZ][ZZ] = z1;
 708         if (sscanf (line, "%*f%*f%*f%lf%lf%lf%lf%lf%lf",
 709                     &x1, &y1, &z1, &x2, &y2, &z2) != 6)
 710         {
 711             x1 = y1 = z1 = x2 = y2 = z2 = 0.0;
 712         }
 713         box[XX][YY] = x1;
 714         box[XX][ZZ] = y1;
 715         box[YY][XX] = z1;
 716         box[YY][ZZ] = x2;
 717         box[ZZ][XX] = y2;
 718         box[ZZ][YY] = z2;
 719     }
 720
 721     return bVel;
 722 }
 723
 724 static void read_whole_conf(const char *infile, char *title,
 725                             t_atoms *atoms, rvec x[], rvec *v, matrix box)
 726 {
 727     FILE    *in;
 728     int      ndec;
 729     t_symtab symtab;
 730
 731     /* open file */
 732     in = gmx_fio_fopen(infile, "r");
 733
 734     open_symtab(&symtab);
 735     get_w_conf(in, infile, title, &symtab, atoms, &ndec, x, v, box);
 736     /* We can't free the symbols, as they are still used in atoms, so
 737      * the only choice is to leak them. */
 738     free_symtab(&symtab);
 739
 740     gmx_fio_fclose(in);
 741 }
 742
 743 static gmx_bool gmx_one_before_eof(FILE *fp)
 744 {
 745     char     data[4];
 746     gmx_bool beof;
 747
 748     if ((beof = fread(data, 1, 1, fp)) == 1)
 749     {
 750         gmx_fseek(fp, -1, SEEK_CUR);
 751     }
 752     return !beof;
 753 }
 754
 755 gmx_bool gro_next_x_or_v(FILE *status, t_trxframe *fr)
 756 {
 757     t_atoms  atoms;
 758     t_symtab symtab;
 759     char     title[STRLEN], *p;
 760     double   tt;
 761     int      ndec = 0, i;
 762
 763     if (gmx_one_before_eof(status))
 764     {
 765         return FALSE;
 766     }
 767
 768     open_symtab(&symtab);
 769     atoms.nr = fr->natoms;
 770     snew(atoms.atom, fr->natoms);
 771     atoms.nres = fr->natoms;
 772     snew(atoms.resinfo, fr->natoms);
 773     snew(atoms.atomname, fr->natoms);
 774
 775     fr->bV    = get_w_conf(status, title, title, &symtab, &atoms, &ndec, fr->x, fr->v, fr->box);
 776     fr->bPrec = TRUE;
 777     fr->prec  = 1;
 778     /* prec = 10^ndec: */
 779     for (i = 0; i < ndec; i++)
 780     {
 781         fr->prec *= 10;
 782     }
 783     fr->title  = title;
 784     fr->bTitle = TRUE;
 785     fr->bX     = TRUE;
 786     fr->bBox   = TRUE;
 787
 788     sfree(atoms.atom);
 789     sfree(atoms.resinfo);
 790     sfree(atoms.atomname);
 791     done_symtab(&symtab);
 792
 793     if ((p = strstr(title, "t=")) != NULL)
 794     {
 795         p += 2;
 796         if (sscanf(p, "%lf", &tt) == 1)
 797         {
 798             fr->time  = tt;
 799             fr->bTime = TRUE;
 800         }
 801         else
 802         {
 803             fr->time  = 0;
 804             fr->bTime = FALSE;
 805         }
 806     }
 807
 808     if (atoms.nr != fr->natoms)
 809     {
 810         gmx_fatal(FARGS, "Number of atoms in gro frame (%d) doesn't match the number in the previous frame (%d)", atoms.nr, fr->natoms);
 811     }
 812
 813     return TRUE;
 814 }
 815
 816 int gro_first_x_or_v(FILE *status, t_trxframe *fr)
 817 {
 818     char title[STRLEN];
 819
 820     frewind(status);
 821     fprintf(stderr, "Reading frames from gro file");
 822     get_coordnum_fp(status, title, &fr->natoms);
 823     frewind(status);
 824     fprintf(stderr, " '%s', %d atoms.\n", title, fr->natoms);
 825     fr->bTitle = TRUE;
 826     fr->title  = title;
 827     if (fr->natoms == 0)
 828     {
 829         gmx_file("No coordinates in gro file");
 830     }
 831
 832     snew(fr->x, fr->natoms);
 833     snew(fr->v, fr->natoms);
 834     gro_next_x_or_v(status, fr);
 835
 836     return fr->natoms;
 837 }
 838
 839 static void make_hconf_format(int pr, gmx_bool bVel, char format[])
 840 {
 841     int l, vpr;
 842
 843     /* build format string for printing,
 844        something like "%8.3f" for x and "%8.4f" for v */
 845     if (pr < 0)
 846     {
 847         pr = 0;
 848     }
 849     if (pr > 30)
 850     {
 851         pr = 30;
 852     }
 853     l   = pr+5;
 854     vpr = pr+1;
 855     if (bVel)
 856     {
 857         sprintf(format, "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
 858                 l, pr, l, pr, l, pr, l, vpr, l, vpr, l, vpr);
 859     }
 860     else
 861     {
 862         sprintf(format, "%%%d.%df%%%d.%df%%%d.%df\n", l, pr, l, pr, l, pr);
 863     }
 864
 865 }
 866
 867 static void write_hconf_box(FILE *out, int pr, matrix box)
 868 {
 869     char format[100];
 870     int  l;
 871
 872     if (pr < 5)
 873     {
 874         pr = 5;
 875     }
 876     l = pr+5;
 877
 878     if (box[XX][YY] || box[XX][ZZ] || box[YY][XX] || box[YY][ZZ] ||
 879         box[ZZ][XX] || box[ZZ][YY])
 880     {
 881         sprintf(format, "%%%d.%df%%%d.%df%%%d.%df"
 882                 "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
 883                 l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr);
 884         fprintf(out, format,
 885                 box[XX][XX], box[YY][YY], box[ZZ][ZZ],
 886                 box[XX][YY], box[XX][ZZ], box[YY][XX],
 887                 box[YY][ZZ], box[ZZ][XX], box[ZZ][YY]);
 888     }
 889     else
 890     {
 891         sprintf(format, "%%%d.%df%%%d.%df%%%d.%df\n", l, pr, l, pr, l, pr);
 892         fprintf(out, format,
 893                 box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
 894     }
 895 }
 896
 897 void write_hconf_indexed_p(FILE *out, const char *title, t_atoms *atoms,
 898                            int nx, const atom_id index[], int pr,
 899                            rvec *x, rvec *v, matrix box)
 900 {
 901     char resnm[6], nm[6], format[100];
 902     int  ai, i, resind, resnr;
 903
 904     bromacs(format, 99);
 905     fprintf (out, "%s\n", (title && title[0]) ? title : format);
 906     fprintf (out, "%5d\n", nx);
 907
 908     make_hconf_format(pr, v != NULL, format);
 909
 910     for (i = 0; (i < nx); i++)
 911     {
 912         ai = index[i];
 913
 914         resind = atoms->atom[ai].resind;
 915         std::strncpy(resnm, " ??? ", sizeof(resnm)-1);
 916         if (resind < atoms->nres)
 917         {
 918             std::strncpy(resnm, *atoms->resinfo[resind].name, sizeof(resnm)-1);
 919             resnr = atoms->resinfo[resind].nr;
 920         }
 921         else
 922         {
 923             std::strncpy(resnm, " ??? ", sizeof(resnm)-1);
 924             resnr = resind + 1;
 925         }
 926
 927         if (atoms->atom)
 928         {
 929             std::strncpy(nm, *atoms->atomname[ai], sizeof(nm)-1);
 930         }
 931         else
 932         {
 933             std::strncpy(nm, " ??? ", sizeof(nm)-1);
 934         }
 935
 936         fprintf(out, "%5d%-5.5s%5.5s%5d", resnr%100000, resnm, nm, (ai+1)%100000);
 937         /* next fprintf uses built format string */
 938         if (v)
 939         {
 940             fprintf(out, format,
 941                     x[ai][XX], x[ai][YY], x[ai][ZZ], v[ai][XX], v[ai][YY], v[ai][ZZ]);
 942         }
 943         else
 944         {
 945             fprintf(out, format,
 946                     x[ai][XX], x[ai][YY], x[ai][ZZ]);
 947         }
 948     }
 949
 950     write_hconf_box(out, pr, box);
 951
 952     fflush(out);
 953 }
 954
 955 static void write_hconf_mtop(FILE *out, const char *title, gmx_mtop_t *mtop,
 956                              int pr,
 957                              rvec *x, rvec *v, matrix box)
 958 {
 959     char                    format[100];
 960     int                     i, resnr;
 961     gmx_mtop_atomloop_all_t aloop;
 962     t_atom                 *atom;
 963     char                   *atomname, *resname;
 964
 965     bromacs(format, 99);
 966     fprintf (out, "%s\n", (title && title[0]) ? title : format);
 967     fprintf (out, "%5d\n", mtop->natoms);
 968
 969     make_hconf_format(pr, v != NULL, format);
 970
 971     aloop = gmx_mtop_atomloop_all_init(mtop);
 972     while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
 973     {
 974         gmx_mtop_atomloop_all_names(aloop, &atomname, &resnr, &resname);
 975
 976         fprintf(out, "%5d%-5.5s%5.5s%5d",
 977                 resnr%100000, resname, atomname, (i+1)%100000);
 978         /* next fprintf uses built format string */
 979         if (v)
 980         {
 981             fprintf(out, format,
 982                     x[i][XX], x[i][YY], x[i][ZZ], v[i][XX], v[i][YY], v[i][ZZ]);
 983         }
 984         else
 985         {
 986             fprintf(out, format,
 987                     x[i][XX], x[i][YY], x[i][ZZ]);
 988         }
 989     }
 990
 991     write_hconf_box(out, pr, box);
 992
 993     fflush(out);
 994 }
 995
 996 void write_hconf_p(FILE *out, const char *title, t_atoms *atoms, int pr,
 997                    rvec *x, rvec *v, matrix box)
 998 {
 999     atom_id *aa;
1000     int      i;
1001
1002     snew(aa, atoms->nr);
1003     for (i = 0; (i < atoms->nr); i++)
1004     {
1005         aa[i] = i;
1006     }
1007     write_hconf_indexed_p(out, title, atoms, atoms->nr, aa, pr, x, v, box);
1008     sfree(aa);
1009 }
1010
1011 void write_conf_p(const char *outfile, const char *title,
1012                   t_atoms *atoms, int pr,
1013                   rvec *x, rvec *v, matrix box)
1014 {
1015     FILE *out;
1016
1017     out = gmx_fio_fopen(outfile, "w");
1018     write_hconf_p(out, title, atoms, pr, x, v, box);
1019
1020     gmx_fio_fclose (out);
1021 }
1022
1023 static void write_conf(const char *outfile, const char *title, t_atoms *atoms,
1024                        rvec *x, rvec *v, matrix box)
1025 {
1026     write_conf_p(outfile, title, atoms, 3, x, v, box);
1027 }
1028
1029 void write_sto_conf_indexed(const char *outfile, const char *title,
1030                             t_atoms *atoms,
1031                             rvec x[], rvec *v, int ePBC, matrix box,
1032                             atom_id nindex, atom_id index[])
1033 {
1034     FILE       *out;
1035     int         ftp;
1036     t_trxframe  fr;
1037
1038     ftp = fn2ftp(outfile);
1039     switch (ftp)
1040     {
1041         case efGRO:
1042             out = gmx_fio_fopen(outfile, "w");
1043             write_hconf_indexed_p(out, title, atoms, nindex, index, 3, x, v, box);
1044             gmx_fio_fclose(out);
1045             break;
1046         case efG96:
1047             clear_trxframe(&fr, TRUE);
1048             fr.bTitle = TRUE;
1049             fr.title  = title;
1050             fr.natoms = atoms->nr;
1051             fr.bAtoms = TRUE;
1052             fr.atoms  = atoms;
1053             fr.bX     = TRUE;
1054             fr.x      = x;
1055             if (v)
1056             {
1057                 fr.bV = TRUE;
1058                 fr.v  = v;
1059             }
1060             fr.bBox = TRUE;
1061             copy_mat(box, fr.box);
1062             out = gmx_fio_fopen(outfile, "w");
1063             write_g96_conf(out, &fr, nindex, index);
1064             gmx_fio_fclose(out);
1065             break;
1066         case efPDB:
1067         case efBRK:
1068         case efENT:
1069         case efPQR:
1070             out = gmx_fio_fopen(outfile, "w");
1071             write_pdbfile_indexed(out, title, atoms, x, ePBC, box, ' ', -1, nindex, index, NULL, TRUE);
1072             gmx_fio_fclose(out);
1073             break;
1074         case efESP:
1075             out = gmx_fio_fopen(outfile, "w");
1076             write_espresso_conf_indexed(out, title, atoms, nindex, index, x, v, box);
1077             gmx_fio_fclose(out);
1078             break;
1079         case efTPR:
1080             gmx_fatal(FARGS, "Sorry, can not write a topology to %s", outfile);
1081             break;
1082         default:
1083             gmx_incons("Not supported in write_sto_conf_indexed");
1084     }
1085 }
1086
1087 void write_sto_conf(const char *outfile, const char *title, t_atoms *atoms,
1088                     rvec x[], rvec *v, int ePBC, matrix box)
1089 {
1090     FILE       *out;
1091     int         ftp;
1092     t_trxframe  fr;
1093
1094     ftp = fn2ftp(outfile);
1095     switch (ftp)
1096     {
1097         case efGRO:
1098             write_conf(outfile, title, atoms, x, v, box);
1099             break;
1100         case efG96:
1101             clear_trxframe(&fr, TRUE);
1102             fr.bTitle = TRUE;
1103             fr.title  = title;
1104             fr.natoms = atoms->nr;
1105             fr.bAtoms = TRUE;
1106             fr.atoms  = atoms;
1107             fr.bX     = TRUE;
1108             fr.x      = x;
1109             if (v)
1110             {
1111                 fr.bV = TRUE;
1112                 fr.v  = v;
1113             }
1114             fr.bBox = TRUE;
1115             copy_mat(box, fr.box);
1116             out = gmx_fio_fopen(outfile, "w");
1117             write_g96_conf(out, &fr, -1, NULL);
1118             gmx_fio_fclose(out);
1119             break;
1120         case efPDB:
1121         case efBRK:
1122         case efENT:
1123             out = gmx_fio_fopen(outfile, "w");
1124             write_pdbfile(out, title, atoms, x, ePBC, box, ' ', -1, NULL, TRUE);
1125             gmx_fio_fclose(out);
1126             break;
1127         case efESP:
1128             out = gmx_fio_fopen(outfile, "w");
1129             write_espresso_conf_indexed(out, title, atoms, atoms->nr, NULL, x, v, box);
1130             gmx_fio_fclose(out);
1131             break;
1132         case efTPR:
1133             gmx_fatal(FARGS, "Sorry, can not write a topology to %s", outfile);
1134             break;
1135         default:
1136             gmx_incons("Not supported in write_sto_conf");
1137     }
1138 }
1139
1140 void write_sto_conf_mtop(const char *outfile, const char *title,
1141                          gmx_mtop_t *mtop,
1142                          rvec x[], rvec *v, int ePBC, matrix box)
1143 {
1144     int     ftp;
1145     FILE   *out;
1146     t_atoms atoms;
1147
1148     ftp = fn2ftp(outfile);
1149     switch (ftp)
1150     {
1151         case efGRO:
1152             out = gmx_fio_fopen(outfile, "w");
1153             write_hconf_mtop(out, title, mtop, 3, x, v, box);
1154             gmx_fio_fclose(out);
1155             break;
1156         default:
1157             /* This is a brute force approach which requires a lot of memory.
1158              * We should implement mtop versions of all writing routines.
1159              */
1160             atoms = gmx_mtop_global_atoms(mtop);
1161
1162             write_sto_conf(outfile, title, &atoms, x, v, ePBC, box);
1163
1164             done_atom(&atoms);
1165             break;
1166     }
1167 }
1168
1169 void get_stx_coordnum(const char *infile, int *natoms)
1170 {
1171     FILE      *in;
1172     int        ftp, tpxver, tpxgen;
1173     t_trxframe fr;
1174     char       g96_line[STRLEN+1];
1175
1176     ftp = fn2ftp(infile);
1177     range_check(ftp, 0, efNR);
1178     switch (ftp)
1179     {
1180         case efGRO:
1181             get_coordnum(infile, natoms);
1182             break;
1183         case efG96:
1184             in        = gmx_fio_fopen(infile, "r");
1185             fr.title  = NULL;
1186             fr.natoms = -1;
1187             fr.atoms  = NULL;
1188             fr.x      = NULL;
1189             fr.v      = NULL;
1190             fr.f      = NULL;
1191             *natoms   = read_g96_conf(in, infile, &fr, g96_line);
1192             gmx_fio_fclose(in);
1193             break;
1194         case efPDB:
1195         case efBRK:
1196         case efENT:
1197             in = gmx_fio_fopen(infile, "r");
1198             get_pdb_coordnum(in, natoms);
1199             gmx_fio_fclose(in);
1200             break;
1201         case efESP:
1202             *natoms = get_espresso_coordnum(infile);
1203             break;
1204         case efTPR:
1205         {
1206             t_tpxheader tpx;
1207
1208             read_tpxheader(infile, &tpx, TRUE, &tpxver, &tpxgen);
1209             *natoms = tpx.natoms;
1210             break;
1211         }
1212         default:
1213             gmx_fatal(FARGS, "File type %s not supported in get_stx_coordnum",
1214                       ftp2ext(ftp));
1215     }
1216 }
1217
1218 static void tpx_make_chain_identifiers(t_atoms *atoms, t_block *mols)
1219 {
1220     int  m, a, a0, a1, r;
1221     char c, chainid;
1222     int  chainnum;
1223
1224     /* We always assign a new chain number, but save the chain id characters
1225      * for larger molecules.
1226      */
1227 #define CHAIN_MIN_ATOMS 15
1228
1229     chainnum = 0;
1230     chainid  = 'A';
1231     for (m = 0; m < mols->nr; m++)
1232     {
1233         a0 = mols->index[m];
1234         a1 = mols->index[m+1];
1235         if ((a1-a0 >= CHAIN_MIN_ATOMS) && (chainid <= 'Z'))
1236         {
1237             c = chainid;
1238             chainid++;
1239         }
1240         else
1241         {
1242             c = ' ';
1243         }
1244         for (a = a0; a < a1; a++)
1245         {
1246             atoms->resinfo[atoms->atom[a].resind].chainnum = chainnum;
1247             atoms->resinfo[atoms->atom[a].resind].chainid  = c;
1248         }
1249         chainnum++;
1250     }
1251
1252     /* Blank out the chain id if there was only one chain */
1253     if (chainid == 'B')
1254     {
1255         for (r = 0; r < atoms->nres; r++)
1256         {
1257             atoms->resinfo[r].chainid = ' ';
1258         }
1259     }
1260 }
1261
1262 void read_stx_conf(const char *infile, char *title, t_atoms *atoms,
1263                    rvec x[], rvec *v, int *ePBC, matrix box)
1264 {
1265     FILE       *in;
1266     gmx_mtop_t *mtop;
1267     t_topology  top;
1268     t_trxframe  fr;
1269     int         i, ftp, natoms;
1270     char        g96_line[STRLEN+1];
1271
1272     if (atoms->nr == 0)
1273     {
1274         fprintf(stderr, "Warning: Number of atoms in %s is 0\n", infile);
1275     }
1276     else if (atoms->atom == NULL)
1277     {
1278         gmx_mem("Uninitialized array atom");
1279     }
1280
1281     if (ePBC)
1282     {
1283         *ePBC = -1;
1284     }
1285
1286     ftp = fn2ftp(infile);
1287     switch (ftp)
1288     {
1289         case efGRO:
1290             read_whole_conf(infile, title, atoms, x, v, box);
1291             break;
1292         case efG96:
1293             fr.title  = NULL;
1294             fr.natoms = atoms->nr;
1295             fr.atoms  = atoms;
1296             fr.x      = x;
1297             fr.v      = v;
1298             fr.f      = NULL;
1299             in        = gmx_fio_fopen(infile, "r");
1300             read_g96_conf(in, infile, &fr, g96_line);
1301             gmx_fio_fclose(in);
1302             copy_mat(fr.box, box);
1303             std::strncpy(title, fr.title, STRLEN);
1304             break;
1305         case efPDB:
1306         case efBRK:
1307         case efENT:
1308             read_pdb_conf(infile, title, atoms, x, ePBC, box, TRUE, NULL);
1309             break;
1310         case efESP:
1311             read_espresso_conf(infile, title, atoms, x, v, box);
1312             break;
1313         case efTPR:
1314             snew(mtop, 1);
1315             i = read_tpx(infile, NULL, box, &natoms, x, v, NULL, mtop);
1316             if (ePBC)
1317             {
1318                 *ePBC = i;
1319             }
1320
1321             strcpy(title, *(mtop->name));
1322
1323             /* Free possibly allocated memory */
1324             done_atom(atoms);
1325
1326             *atoms = gmx_mtop_global_atoms(mtop);
1327             top    = gmx_mtop_t_to_t_topology(mtop);
1328             tpx_make_chain_identifiers(atoms, &top.mols);
1329
1330             sfree(mtop);
1331             /* The strings in the symtab are still in use in the returned t_atoms
1332              * structure, so we should not free them. But there is no place to put the
1333              * symbols; the only choice is to leak the memory...
1334              * So we clear the symbol table before freeing the topology structure. */
1335             free_symtab(&top.symtab);
1336             done_top(&top);
1337
1338             break;
1339         default:
1340             gmx_incons("Not supported in read_stx_conf");
1341     }
1342 }
1343
1344 static void done_gmx_groups_t(gmx_groups_t *g)
1345 {
1346     int i;
1347
1348     for (i = 0; (i < egcNR); i++)
1349     {
1350         if (NULL != g->grps[i].nm_ind)
1351         {
1352             sfree(g->grps[i].nm_ind);
1353             g->grps[i].nm_ind = NULL;
1354         }
1355         if (NULL != g->grpnr[i])
1356         {
1357             sfree(g->grpnr[i]);
1358             g->grpnr[i] = NULL;
1359         }
1360     }
1361     /* The contents of this array is in symtab, don't free it here */
1362     sfree(g->grpname);
1363 }
1364
1365 gmx_bool read_tps_conf(const char *infile, char *title, t_topology *top, int *ePBC,
1366                        rvec **x, rvec **v, matrix box, gmx_bool bMass)
1367 {
1368     t_tpxheader      header;
1369     int              natoms, i, version, generation;
1370     gmx_bool         bTop, bXNULL = FALSE;
1371     gmx_mtop_t      *mtop;
1372     gmx_atomprop_t   aps;
1373
1374     bTop  = fn2bTPX(infile);
1375     *ePBC = -1;
1376     if (bTop)
1377     {
1378         read_tpxheader(infile, &header, TRUE, &version, &generation);
1379         if (x)
1380         {
1381             snew(*x, header.natoms);
1382         }
1383         if (v)
1384         {
1385             snew(*v, header.natoms);
1386         }
1387         snew(mtop, 1);
1388         *ePBC = read_tpx(infile, NULL, box, &natoms,
1389                          (x == NULL) ? NULL : *x, (v == NULL) ? NULL : *v, NULL, mtop);
1390         *top = gmx_mtop_t_to_t_topology(mtop);
1391         /* In this case we need to throw away the group data too */
1392         done_gmx_groups_t(&mtop->groups);
1393         sfree(mtop);
1394         std::strcpy(title, *top->name);
1395         tpx_make_chain_identifiers(&top->atoms, &top->mols);
1396     }
1397     else
1398     {
1399         get_stx_coordnum(infile, &natoms);
1400         init_t_atoms(&top->atoms, natoms, (fn2ftp(infile) == efPDB));
1401         if (x == NULL)
1402         {
1403             snew(x, 1);
1404             bXNULL = TRUE;
1405         }
1406         snew(*x, natoms);
1407         if (v)
1408         {
1409             snew(*v, natoms);
1410         }
1411         read_stx_conf(infile, title, &top->atoms, *x, (v == NULL) ? NULL : *v, ePBC, box);
1412         if (bXNULL)
1413         {
1414             sfree(*x);
1415             sfree(x);
1416         }
1417         if (bMass)
1418         {
1419             aps = gmx_atomprop_init();
1420             for (i = 0; (i < natoms); i++)
1421             {
1422                 if (!gmx_atomprop_query(aps, epropMass,
1423                                         *top->atoms.resinfo[top->atoms.atom[i].resind].name,
1424                                         *top->atoms.atomname[i],
1425                                         &(top->atoms.atom[i].m)))
1426                 {
1427                     if (debug)
1428                     {
1429                         fprintf(debug, "Can not find mass for atom %s %d %s, setting to 1\n",
1430                                 *top->atoms.resinfo[top->atoms.atom[i].resind].name,
1431                                 top->atoms.resinfo[top->atoms.atom[i].resind].nr,
1432                                 *top->atoms.atomname[i]);
1433                     }
1434                 }
1435             }
1436             gmx_atomprop_destroy(aps);
1437         }
1438         top->idef.ntypes = -1;
1439     }
1440
1441     return bTop;
1442 }