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/groio.h"
  53 #include "gromacs/fileio/pdbio.h"
  54 #include "gromacs/fileio/tpxio.h"
  55 #include "gromacs/fileio/trx.h"
  56 #include "gromacs/fileio/trxio.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 void write_sto_conf_indexed(const char *outfile, const char *title,
 487                             t_atoms *atoms,
 488                             rvec x[], rvec *v, int ePBC, matrix box,
 489                             atom_id nindex, atom_id index[])
 490 {
 491     FILE       *out;
 492     int         ftp;
 493     t_trxframe  fr;
 494
 495     ftp = fn2ftp(outfile);
 496     switch (ftp)
 497     {
 498         case efGRO:
 499             out = gmx_fio_fopen(outfile, "w");
 500             write_hconf_indexed_p(out, title, atoms, nindex, index, 3, x, v, box);
 501             gmx_fio_fclose(out);
 502             break;
 503         case efG96:
 504             clear_trxframe(&fr, TRUE);
 505             fr.bTitle = TRUE;
 506             fr.title  = title;
 507             fr.natoms = atoms->nr;
 508             fr.bAtoms = TRUE;
 509             fr.atoms  = atoms;
 510             fr.bX     = TRUE;
 511             fr.x      = x;
 512             if (v)
 513             {
 514                 fr.bV = TRUE;
 515                 fr.v  = v;
 516             }
 517             fr.bBox = TRUE;
 518             copy_mat(box, fr.box);
 519             out = gmx_fio_fopen(outfile, "w");
 520             write_g96_conf(out, &fr, nindex, index);
 521             gmx_fio_fclose(out);
 522             break;
 523         case efPDB:
 524         case efBRK:
 525         case efENT:
 526         case efPQR:
 527             out = gmx_fio_fopen(outfile, "w");
 528             write_pdbfile_indexed(out, title, atoms, x, ePBC, box, ' ', -1, nindex, index, NULL, TRUE);
 529             gmx_fio_fclose(out);
 530             break;
 531         case efESP:
 532             out = gmx_fio_fopen(outfile, "w");
 533             write_espresso_conf_indexed(out, title, atoms, nindex, index, x, v, box);
 534             gmx_fio_fclose(out);
 535             break;
 536         case efTPR:
 537             gmx_fatal(FARGS, "Sorry, can not write a topology to %s", outfile);
 538             break;
 539         default:
 540             gmx_incons("Not supported in write_sto_conf_indexed");
 541     }
 542 }
 543
 544 void write_sto_conf(const char *outfile, const char *title, t_atoms *atoms,
 545                     rvec x[], rvec *v, int ePBC, matrix box)
 546 {
 547     FILE       *out;
 548     int         ftp;
 549     t_trxframe  fr;
 550
 551     ftp = fn2ftp(outfile);
 552     switch (ftp)
 553     {
 554         case efGRO:
 555             write_conf_p(outfile, title, atoms, 3, x, v, box);
 556             break;
 557         case efG96:
 558             clear_trxframe(&fr, TRUE);
 559             fr.bTitle = TRUE;
 560             fr.title  = title;
 561             fr.natoms = atoms->nr;
 562             fr.bAtoms = TRUE;
 563             fr.atoms  = atoms;
 564             fr.bX     = TRUE;
 565             fr.x      = x;
 566             if (v)
 567             {
 568                 fr.bV = TRUE;
 569                 fr.v  = v;
 570             }
 571             fr.bBox = TRUE;
 572             copy_mat(box, fr.box);
 573             out = gmx_fio_fopen(outfile, "w");
 574             write_g96_conf(out, &fr, -1, NULL);
 575             gmx_fio_fclose(out);
 576             break;
 577         case efPDB:
 578         case efBRK:
 579         case efENT:
 580             out = gmx_fio_fopen(outfile, "w");
 581             write_pdbfile(out, title, atoms, x, ePBC, box, ' ', -1, NULL, TRUE);
 582             gmx_fio_fclose(out);
 583             break;
 584         case efESP:
 585             out = gmx_fio_fopen(outfile, "w");
 586             write_espresso_conf_indexed(out, title, atoms, atoms->nr, NULL, x, v, box);
 587             gmx_fio_fclose(out);
 588             break;
 589         case efTPR:
 590             gmx_fatal(FARGS, "Sorry, can not write a topology to %s", outfile);
 591             break;
 592         default:
 593             gmx_incons("Not supported in write_sto_conf");
 594     }
 595 }
 596
 597 void write_sto_conf_mtop(const char *outfile, const char *title,
 598                          gmx_mtop_t *mtop,
 599                          rvec x[], rvec *v, int ePBC, matrix box)
 600 {
 601     int     ftp;
 602     FILE   *out;
 603     t_atoms atoms;
 604
 605     ftp = fn2ftp(outfile);
 606     switch (ftp)
 607     {
 608         case efGRO:
 609             out = gmx_fio_fopen(outfile, "w");
 610             write_hconf_mtop(out, title, mtop, 3, x, v, box);
 611             gmx_fio_fclose(out);
 612             break;
 613         default:
 614             /* This is a brute force approach which requires a lot of memory.
 615              * We should implement mtop versions of all writing routines.
 616              */
 617             atoms = gmx_mtop_global_atoms(mtop);
 618
 619             write_sto_conf(outfile, title, &atoms, x, v, ePBC, box);
 620
 621             done_atom(&atoms);
 622             break;
 623     }
 624 }
 625
 626 void get_stx_coordnum(const char *infile, int *natoms)
 627 {
 628     FILE      *in;
 629     int        ftp, tpxver, tpxgen;
 630     t_trxframe fr;
 631     char       g96_line[STRLEN+1];
 632
 633     ftp = fn2ftp(infile);
 634     range_check(ftp, 0, efNR);
 635     switch (ftp)
 636     {
 637         case efGRO:
 638             get_coordnum(infile, natoms);
 639             break;
 640         case efG96:
 641             in        = gmx_fio_fopen(infile, "r");
 642             fr.title  = NULL;
 643             fr.natoms = -1;
 644             fr.atoms  = NULL;
 645             fr.x      = NULL;
 646             fr.v      = NULL;
 647             fr.f      = NULL;
 648             *natoms   = read_g96_conf(in, infile, &fr, g96_line);
 649             gmx_fio_fclose(in);
 650             break;
 651         case efPDB:
 652         case efBRK:
 653         case efENT:
 654             in = gmx_fio_fopen(infile, "r");
 655             get_pdb_coordnum(in, natoms);
 656             gmx_fio_fclose(in);
 657             break;
 658         case efESP:
 659             *natoms = get_espresso_coordnum(infile);
 660             break;
 661         case efTPR:
 662         {
 663             t_tpxheader tpx;
 664
 665             read_tpxheader(infile, &tpx, TRUE, &tpxver, &tpxgen);
 666             *natoms = tpx.natoms;
 667             break;
 668         }
 669         default:
 670             gmx_fatal(FARGS, "File type %s not supported in get_stx_coordnum",
 671                       ftp2ext(ftp));
 672     }
 673 }
 674
 675 static void tpx_make_chain_identifiers(t_atoms *atoms, t_block *mols)
 676 {
 677     int  m, a, a0, a1, r;
 678     char c, chainid;
 679     int  chainnum;
 680
 681     /* We always assign a new chain number, but save the chain id characters
 682      * for larger molecules.
 683      */
 684 #define CHAIN_MIN_ATOMS 15
 685
 686     chainnum = 0;
 687     chainid  = 'A';
 688     for (m = 0; m < mols->nr; m++)
 689     {
 690         a0 = mols->index[m];
 691         a1 = mols->index[m+1];
 692         if ((a1-a0 >= CHAIN_MIN_ATOMS) && (chainid <= 'Z'))
 693         {
 694             c = chainid;
 695             chainid++;
 696         }
 697         else
 698         {
 699             c = ' ';
 700         }
 701         for (a = a0; a < a1; a++)
 702         {
 703             atoms->resinfo[atoms->atom[a].resind].chainnum = chainnum;
 704             atoms->resinfo[atoms->atom[a].resind].chainid  = c;
 705         }
 706         chainnum++;
 707     }
 708
 709     /* Blank out the chain id if there was only one chain */
 710     if (chainid == 'B')
 711     {
 712         for (r = 0; r < atoms->nres; r++)
 713         {
 714             atoms->resinfo[r].chainid = ' ';
 715         }
 716     }
 717 }
 718
 719 void read_stx_conf(const char *infile, char *title, t_atoms *atoms,
 720                    rvec x[], rvec *v, int *ePBC, matrix box)
 721 {
 722     FILE       *in;
 723     gmx_mtop_t *mtop;
 724     t_topology  top;
 725     t_trxframe  fr;
 726     int         i, ftp, natoms;
 727     char        g96_line[STRLEN+1];
 728
 729     if (atoms->nr == 0)
 730     {
 731         fprintf(stderr, "Warning: Number of atoms in %s is 0\n", infile);
 732     }
 733     else if (atoms->atom == NULL)
 734     {
 735         gmx_mem("Uninitialized array atom");
 736     }
 737
 738     if (ePBC)
 739     {
 740         *ePBC = -1;
 741     }
 742
 743     ftp = fn2ftp(infile);
 744     switch (ftp)
 745     {
 746         case efGRO:
 747             read_whole_conf(infile, title, atoms, x, v, box);
 748             break;
 749         case efG96:
 750             fr.title  = NULL;
 751             fr.natoms = atoms->nr;
 752             fr.atoms  = atoms;
 753             fr.x      = x;
 754             fr.v      = v;
 755             fr.f      = NULL;
 756             in        = gmx_fio_fopen(infile, "r");
 757             read_g96_conf(in, infile, &fr, g96_line);
 758             gmx_fio_fclose(in);
 759             copy_mat(fr.box, box);
 760             std::strncpy(title, fr.title, STRLEN);
 761             break;
 762         case efPDB:
 763         case efBRK:
 764         case efENT:
 765             read_pdb_conf(infile, title, atoms, x, ePBC, box, TRUE, NULL);
 766             break;
 767         case efESP:
 768             read_espresso_conf(infile, title, atoms, x, v, box);
 769             break;
 770         case efTPR:
 771             snew(mtop, 1);
 772             i = read_tpx(infile, NULL, box, &natoms, x, v, NULL, mtop);
 773             if (ePBC)
 774             {
 775                 *ePBC = i;
 776             }
 777
 778             strcpy(title, *(mtop->name));
 779
 780             /* Free possibly allocated memory */
 781             done_atom(atoms);
 782
 783             *atoms = gmx_mtop_global_atoms(mtop);
 784             top    = gmx_mtop_t_to_t_topology(mtop);
 785             tpx_make_chain_identifiers(atoms, &top.mols);
 786
 787             sfree(mtop);
 788             /* The strings in the symtab are still in use in the returned t_atoms
 789              * structure, so we should not free them. But there is no place to put the
 790              * symbols; the only choice is to leak the memory...
 791              * So we clear the symbol table before freeing the topology structure. */
 792             free_symtab(&top.symtab);
 793             done_top(&top);
 794
 795             break;
 796         default:
 797             gmx_incons("Not supported in read_stx_conf");
 798     }
 799 }
 800
 801 static void done_gmx_groups_t(gmx_groups_t *g)
 802 {
 803     int i;
 804
 805     for (i = 0; (i < egcNR); i++)
 806     {
 807         if (NULL != g->grps[i].nm_ind)
 808         {
 809             sfree(g->grps[i].nm_ind);
 810             g->grps[i].nm_ind = NULL;
 811         }
 812         if (NULL != g->grpnr[i])
 813         {
 814             sfree(g->grpnr[i]);
 815             g->grpnr[i] = NULL;
 816         }
 817     }
 818     /* The contents of this array is in symtab, don't free it here */
 819     sfree(g->grpname);
 820 }
 821
 822 gmx_bool read_tps_conf(const char *infile, char *title, t_topology *top, int *ePBC,
 823                        rvec **x, rvec **v, matrix box, gmx_bool bMass)
 824 {
 825     t_tpxheader      header;
 826     int              natoms, i, version, generation;
 827     gmx_bool         bTop, bXNULL = FALSE;
 828     gmx_mtop_t      *mtop;
 829     gmx_atomprop_t   aps;
 830
 831     bTop  = fn2bTPX(infile);
 832     *ePBC = -1;
 833     if (bTop)
 834     {
 835         read_tpxheader(infile, &header, TRUE, &version, &generation);
 836         if (x)
 837         {
 838             snew(*x, header.natoms);
 839         }
 840         if (v)
 841         {
 842             snew(*v, header.natoms);
 843         }
 844         snew(mtop, 1);
 845         *ePBC = read_tpx(infile, NULL, box, &natoms,
 846                          (x == NULL) ? NULL : *x, (v == NULL) ? NULL : *v, NULL, mtop);
 847         *top = gmx_mtop_t_to_t_topology(mtop);
 848         /* In this case we need to throw away the group data too */
 849         done_gmx_groups_t(&mtop->groups);
 850         sfree(mtop);
 851         std::strcpy(title, *top->name);
 852         tpx_make_chain_identifiers(&top->atoms, &top->mols);
 853     }
 854     else
 855     {
 856         get_stx_coordnum(infile, &natoms);
 857         init_t_atoms(&top->atoms, natoms, (fn2ftp(infile) == efPDB));
 858         if (x == NULL)
 859         {
 860             snew(x, 1);
 861             bXNULL = TRUE;
 862         }
 863         snew(*x, natoms);
 864         if (v)
 865         {
 866             snew(*v, natoms);
 867         }
 868         read_stx_conf(infile, title, &top->atoms, *x, (v == NULL) ? NULL : *v, ePBC, box);
 869         if (bXNULL)
 870         {
 871             sfree(*x);
 872             sfree(x);
 873         }
 874         if (bMass)
 875         {
 876             aps = gmx_atomprop_init();
 877             for (i = 0; (i < natoms); i++)
 878             {
 879                 if (!gmx_atomprop_query(aps, epropMass,
 880                                         *top->atoms.resinfo[top->atoms.atom[i].resind].name,
 881                                         *top->atoms.atomname[i],
 882                                         &(top->atoms.atom[i].m)))
 883                 {
 884                     if (debug)
 885                     {
 886                         fprintf(debug, "Can not find mass for atom %s %d %s, setting to 1\n",
 887                                 *top->atoms.resinfo[top->atoms.atom[i].resind].name,
 888                                 top->atoms.resinfo[top->atoms.atom[i].resind].nr,
 889                                 *top->atoms.atomname[i]);
 890                     }
 891                 }
 892             }
 893             gmx_atomprop_destroy(aps);
 894         }
 895         top->idef.ntypes = -1;
 896     }
 897
 898     return bTop;
 899 }