src/gromacs/fileio/groio.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 "groio.h"
  40
  41 #include <cstdio>
  42 #include <cstring>
  43
  44 #include <algorithm>
  45
  46 #include "gromacs/fileio/gmxfio.h"
  47 #include "gromacs/fileio/trx.h"
  48 #include "gromacs/legacyheaders/copyrite.h"
  49 #include "gromacs/topology/atoms.h"
  50 #include "gromacs/topology/mtop_util.h"
  51 #include "gromacs/topology/symtab.h"
  52 #include "gromacs/topology/topology.h"
  53 #include "gromacs/utility/cstringutil.h"
  54 #include "gromacs/utility/fatalerror.h"
  55 #include "gromacs/utility/futil.h"
  56 #include "gromacs/utility/smalloc.h"
  57
  58 static void get_coordnum_fp(FILE *in, char *title, int *natoms)
  59 {
  60     char line[STRLEN+1];
  61
  62     fgets2(title, STRLEN, in);
  63     fgets2(line, STRLEN, in);
  64     if (sscanf(line, "%d", natoms) != 1)
  65     {
  66         gmx_fatal(FARGS, "gro file does not have the number of atoms on the second line");
  67     }
  68 }
  69
  70 void get_coordnum(const char *infile, int *natoms)
  71 {
  72     FILE *in;
  73     char  title[STRLEN];
  74
  75     in = gmx_fio_fopen(infile, "r");
  76     get_coordnum_fp(in, title, natoms);
  77     gmx_fio_fclose(in);
  78 }
  79
  80 static gmx_bool get_w_conf(FILE *in, const char *infile, char *title,
  81                            t_symtab *symtab, t_atoms *atoms, int *ndec,
  82                            rvec x[], rvec *v, matrix box)
  83 {
  84     char       name[6];
  85     char       resname[6], oldresname[6];
  86     char       line[STRLEN+1], *ptr;
  87     char       buf[256];
  88     double     x1, y1, z1, x2, y2, z2;
  89     rvec       xmin, xmax;
  90     int        natoms, i, m, resnr, newres, oldres, ddist, c;
  91     gmx_bool   bFirst, bVel;
  92     char      *p1, *p2, *p3;
  93
  94     newres  = -1;
  95     oldres  = -12345; /* Unlikely number for the first residue! */
  96     ddist   = 0;
  97
  98     /* Read the title and number of atoms */
  99     get_coordnum_fp(in, title, &natoms);
 100
 101     if (natoms > atoms->nr)
 102     {
 103         gmx_fatal(FARGS, "gro file contains more atoms (%d) than expected (%d)",
 104                   natoms, atoms->nr);
 105     }
 106     else if (natoms <  atoms->nr)
 107     {
 108         fprintf(stderr, "Warning: gro file contains less atoms (%d) than expected"
 109                 " (%d)\n", natoms, atoms->nr);
 110     }
 111
 112     bFirst = TRUE;
 113
 114     bVel = FALSE;
 115
 116     resname[0]     = '\0';
 117     oldresname[0]  = '\0';
 118
 119     /* just pray the arrays are big enough */
 120     for (i = 0; (i < natoms); i++)
 121     {
 122         if ((fgets2(line, STRLEN, in)) == NULL)
 123         {
 124             gmx_fatal(FARGS, "Unexpected end of file in file %s at line %d",
 125                       infile, i+2);
 126         }
 127         if (strlen(line) < 39)
 128         {
 129             gmx_fatal(FARGS, "Invalid line in %s for atom %d:\n%s", infile, i+1, line);
 130         }
 131
 132         /* determine read precision from distance between periods
 133            (decimal points) */
 134         if (bFirst)
 135         {
 136             bFirst = FALSE;
 137             p1     = strchr(line, '.');
 138             if (p1 == NULL)
 139             {
 140                 gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
 141             }
 142             p2 = strchr(&p1[1], '.');
 143             if (p2 == NULL)
 144             {
 145                 gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
 146             }
 147             ddist = p2 - p1;
 148             *ndec = ddist - 5;
 149
 150             p3 = strchr(&p2[1], '.');
 151             if (p3 == NULL)
 152             {
 153                 gmx_fatal(FARGS, "A coordinate in file %s does not contain a '.'", infile);
 154             }
 155
 156             if (p3 - p2 != ddist)
 157             {
 158                 gmx_fatal(FARGS, "The spacing of the decimal points in file %s is not consistent for x, y and z", infile);
 159             }
 160         }
 161
 162         /* residue number*/
 163         memcpy(name, line, 5);
 164         name[5] = '\0';
 165         sscanf(name, "%d", &resnr);
 166         sscanf(line+5, "%5s", resname);
 167
 168         if (resnr != oldres || strncmp(resname, oldresname, sizeof(resname)))
 169         {
 170             oldres = resnr;
 171             newres++;
 172             if (newres >= natoms)
 173             {
 174                 gmx_fatal(FARGS, "More residues than atoms in %s (natoms = %d)",
 175                           infile, natoms);
 176             }
 177             atoms->atom[i].resind = newres;
 178             t_atoms_set_resinfo(atoms, i, symtab, resname, resnr, ' ', 0, ' ');
 179         }
 180         else
 181         {
 182             atoms->atom[i].resind = newres;
 183         }
 184
 185         /* atomname */
 186         std::memcpy(name, line+10, 5);
 187         atoms->atomname[i] = put_symtab(symtab, name);
 188
 189         /* Copy resname to oldresname after we are done with the sanity check above */
 190         std::strncpy(oldresname, resname, sizeof(oldresname));
 191
 192         /* eventueel controle atomnumber met i+1 */
 193
 194         /* coordinates (start after residue data) */
 195         ptr = line + 20;
 196         /* Read fixed format */
 197         for (m = 0; m < DIM; m++)
 198         {
 199             for (c = 0; (c < ddist && ptr[0]); c++)
 200             {
 201                 buf[c] = ptr[0];
 202                 ptr++;
 203             }
 204             buf[c] = '\0';
 205             if (sscanf(buf, "%lf %lf", &x1, &x2) != 1)
 206             {
 207                 gmx_fatal(FARGS, "Something is wrong in the coordinate formatting of file %s. Note that gro is fixed format (see the manual)", infile);
 208             }
 209             else
 210             {
 211                 x[i][m] = x1;
 212             }
 213         }
 214
 215         /* velocities (start after residues and coordinates) */
 216         if (v)
 217         {
 218             /* Read fixed format */
 219             for (m = 0; m < DIM; m++)
 220             {
 221                 for (c = 0; (c < ddist && ptr[0]); c++)
 222                 {
 223                     buf[c] = ptr[0];
 224                     ptr++;
 225                 }
 226                 buf[c] = '\0';
 227                 if (sscanf(buf, "%lf", &x1) != 1)
 228                 {
 229                     v[i][m] = 0;
 230                 }
 231                 else
 232                 {
 233                     v[i][m] = x1;
 234                     bVel    = TRUE;
 235                 }
 236             }
 237         }
 238     }
 239     atoms->nres = newres + 1;
 240
 241     /* box */
 242     fgets2(line, STRLEN, in);
 243     if (sscanf(line, "%lf%lf%lf", &x1, &y1, &z1) != 3)
 244     {
 245         gmx_warning("Bad box in file %s", infile);
 246
 247         /* Generate a cubic box */
 248         for (m = 0; (m < DIM); m++)
 249         {
 250             xmin[m] = xmax[m] = x[0][m];
 251         }
 252         for (i = 1; (i < atoms->nr); i++)
 253         {
 254             for (m = 0; (m < DIM); m++)
 255             {
 256                 xmin[m] = std::min(xmin[m], x[i][m]);
 257                 xmax[m] = std::max(xmax[m], x[i][m]);
 258             }
 259         }
 260         for (i = 0; i < DIM; i++)
 261         {
 262             for (m = 0; m < DIM; m++)
 263             {
 264                 box[i][m] = 0.0;
 265             }
 266         }
 267         for (m = 0; (m < DIM); m++)
 268         {
 269             box[m][m] = (xmax[m]-xmin[m]);
 270         }
 271         fprintf(stderr, "Generated a cubic box %8.3f x %8.3f x %8.3f\n",
 272                 box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
 273     }
 274     else
 275     {
 276         /* We found the first three values, the diagonal elements */
 277         box[XX][XX] = x1;
 278         box[YY][YY] = y1;
 279         box[ZZ][ZZ] = z1;
 280         if (sscanf (line, "%*f%*f%*f%lf%lf%lf%lf%lf%lf",
 281                     &x1, &y1, &z1, &x2, &y2, &z2) != 6)
 282         {
 283             x1 = y1 = z1 = x2 = y2 = z2 = 0.0;
 284         }
 285         box[XX][YY] = x1;
 286         box[XX][ZZ] = y1;
 287         box[YY][XX] = z1;
 288         box[YY][ZZ] = x2;
 289         box[ZZ][XX] = y2;
 290         box[ZZ][YY] = z2;
 291     }
 292
 293     return bVel;
 294 }
 295
 296 void read_whole_conf(const char *infile, char *title,
 297                      t_atoms *atoms, rvec x[], rvec *v, matrix box)
 298 {
 299     FILE    *in;
 300     int      ndec;
 301     t_symtab symtab;
 302
 303     /* open file */
 304     in = gmx_fio_fopen(infile, "r");
 305
 306     open_symtab(&symtab);
 307     get_w_conf(in, infile, title, &symtab, atoms, &ndec, x, v, box);
 308     /* We can't free the symbols, as they are still used in atoms, so
 309      * the only choice is to leak them. */
 310     free_symtab(&symtab);
 311
 312     gmx_fio_fclose(in);
 313 }
 314
 315 static gmx_bool gmx_one_before_eof(FILE *fp)
 316 {
 317     char     data[4];
 318     gmx_bool beof;
 319
 320     if ((beof = fread(data, 1, 1, fp)) == 1)
 321     {
 322         gmx_fseek(fp, -1, SEEK_CUR);
 323     }
 324     return !beof;
 325 }
 326
 327 gmx_bool gro_next_x_or_v(FILE *status, t_trxframe *fr)
 328 {
 329     t_atoms  atoms;
 330     t_symtab symtab;
 331     char     title[STRLEN], *p;
 332     double   tt;
 333     int      ndec = 0, i;
 334
 335     if (gmx_one_before_eof(status))
 336     {
 337         return FALSE;
 338     }
 339
 340     open_symtab(&symtab);
 341     atoms.nr = fr->natoms;
 342     snew(atoms.atom, fr->natoms);
 343     atoms.nres = fr->natoms;
 344     snew(atoms.resinfo, fr->natoms);
 345     snew(atoms.atomname, fr->natoms);
 346
 347     fr->bV    = get_w_conf(status, title, title, &symtab, &atoms, &ndec, fr->x, fr->v, fr->box);
 348     fr->bPrec = TRUE;
 349     fr->prec  = 1;
 350     /* prec = 10^ndec: */
 351     for (i = 0; i < ndec; i++)
 352     {
 353         fr->prec *= 10;
 354     }
 355     fr->title  = title;
 356     fr->bTitle = TRUE;
 357     fr->bX     = TRUE;
 358     fr->bBox   = TRUE;
 359
 360     sfree(atoms.atom);
 361     sfree(atoms.resinfo);
 362     sfree(atoms.atomname);
 363     done_symtab(&symtab);
 364
 365     if ((p = strstr(title, "t=")) != NULL)
 366     {
 367         p += 2;
 368         if (sscanf(p, "%lf", &tt) == 1)
 369         {
 370             fr->time  = tt;
 371             fr->bTime = TRUE;
 372         }
 373         else
 374         {
 375             fr->time  = 0;
 376             fr->bTime = FALSE;
 377         }
 378     }
 379
 380     if (atoms.nr != fr->natoms)
 381     {
 382         gmx_fatal(FARGS, "Number of atoms in gro frame (%d) doesn't match the number in the previous frame (%d)", atoms.nr, fr->natoms);
 383     }
 384
 385     return TRUE;
 386 }
 387
 388 int gro_first_x_or_v(FILE *status, t_trxframe *fr)
 389 {
 390     char title[STRLEN];
 391
 392     frewind(status);
 393     fprintf(stderr, "Reading frames from gro file");
 394     get_coordnum_fp(status, title, &fr->natoms);
 395     frewind(status);
 396     fprintf(stderr, " '%s', %d atoms.\n", title, fr->natoms);
 397     fr->bTitle = TRUE;
 398     fr->title  = title;
 399     if (fr->natoms == 0)
 400     {
 401         gmx_file("No coordinates in gro file");
 402     }
 403
 404     snew(fr->x, fr->natoms);
 405     snew(fr->v, fr->natoms);
 406     gro_next_x_or_v(status, fr);
 407
 408     return fr->natoms;
 409 }
 410
 411 static void make_hconf_format(int pr, gmx_bool bVel, char format[])
 412 {
 413     int l, vpr;
 414
 415     /* build format string for printing,
 416        something like "%8.3f" for x and "%8.4f" for v */
 417     if (pr < 0)
 418     {
 419         pr = 0;
 420     }
 421     if (pr > 30)
 422     {
 423         pr = 30;
 424     }
 425     l   = pr+5;
 426     vpr = pr+1;
 427     if (bVel)
 428     {
 429         sprintf(format, "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
 430                 l, pr, l, pr, l, pr, l, vpr, l, vpr, l, vpr);
 431     }
 432     else
 433     {
 434         sprintf(format, "%%%d.%df%%%d.%df%%%d.%df\n", l, pr, l, pr, l, pr);
 435     }
 436
 437 }
 438
 439 static void write_hconf_box(FILE *out, int pr, matrix box)
 440 {
 441     char format[100];
 442     int  l;
 443
 444     if (pr < 5)
 445     {
 446         pr = 5;
 447     }
 448     l = pr+5;
 449
 450     if (box[XX][YY] || box[XX][ZZ] || box[YY][XX] || box[YY][ZZ] ||
 451         box[ZZ][XX] || box[ZZ][YY])
 452     {
 453         sprintf(format, "%%%d.%df%%%d.%df%%%d.%df"
 454                 "%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df%%%d.%df\n",
 455                 l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr, l, pr);
 456         fprintf(out, format,
 457                 box[XX][XX], box[YY][YY], box[ZZ][ZZ],
 458                 box[XX][YY], box[XX][ZZ], box[YY][XX],
 459                 box[YY][ZZ], box[ZZ][XX], box[ZZ][YY]);
 460     }
 461     else
 462     {
 463         sprintf(format, "%%%d.%df%%%d.%df%%%d.%df\n", l, pr, l, pr, l, pr);
 464         fprintf(out, format,
 465                 box[XX][XX], box[YY][YY], box[ZZ][ZZ]);
 466     }
 467 }
 468
 469 void write_hconf_indexed_p(FILE *out, const char *title, t_atoms *atoms,
 470                            int nx, const atom_id index[], int pr,
 471                            rvec *x, rvec *v, matrix box)
 472 {
 473     char resnm[6], nm[6], format[100];
 474     int  ai, i, resind, resnr;
 475
 476     bromacs(format, 99);
 477     fprintf(out, "%s\n", (title && title[0]) ? title : format);
 478     fprintf(out, "%5d\n", nx);
 479
 480     make_hconf_format(pr, v != NULL, format);
 481
 482     for (i = 0; (i < nx); i++)
 483     {
 484         ai = index[i];
 485
 486         resind = atoms->atom[ai].resind;
 487         std::strncpy(resnm, " ??? ", sizeof(resnm)-1);
 488         if (resind < atoms->nres)
 489         {
 490             std::strncpy(resnm, *atoms->resinfo[resind].name, sizeof(resnm)-1);
 491             resnr = atoms->resinfo[resind].nr;
 492         }
 493         else
 494         {
 495             std::strncpy(resnm, " ??? ", sizeof(resnm)-1);
 496             resnr = resind + 1;
 497         }
 498
 499         if (atoms->atom)
 500         {
 501             std::strncpy(nm, *atoms->atomname[ai], sizeof(nm)-1);
 502         }
 503         else
 504         {
 505             std::strncpy(nm, " ??? ", sizeof(nm)-1);
 506         }
 507
 508         fprintf(out, "%5d%-5.5s%5.5s%5d", resnr%100000, resnm, nm, (ai+1)%100000);
 509         /* next fprintf uses built format string */
 510         if (v)
 511         {
 512             fprintf(out, format,
 513                     x[ai][XX], x[ai][YY], x[ai][ZZ], v[ai][XX], v[ai][YY], v[ai][ZZ]);
 514         }
 515         else
 516         {
 517             fprintf(out, format,
 518                     x[ai][XX], x[ai][YY], x[ai][ZZ]);
 519         }
 520     }
 521
 522     write_hconf_box(out, pr, box);
 523
 524     fflush(out);
 525 }
 526
 527 void write_hconf_mtop(FILE *out, const char *title, gmx_mtop_t *mtop, int pr,
 528                       rvec *x, rvec *v, matrix box)
 529 {
 530     char                    format[100];
 531     int                     i, resnr;
 532     gmx_mtop_atomloop_all_t aloop;
 533     t_atom                 *atom;
 534     char                   *atomname, *resname;
 535
 536     bromacs(format, 99);
 537     fprintf(out, "%s\n", (title && title[0]) ? title : format);
 538     fprintf(out, "%5d\n", mtop->natoms);
 539
 540     make_hconf_format(pr, v != NULL, format);
 541
 542     aloop = gmx_mtop_atomloop_all_init(mtop);
 543     while (gmx_mtop_atomloop_all_next(aloop, &i, &atom))
 544     {
 545         gmx_mtop_atomloop_all_names(aloop, &atomname, &resnr, &resname);
 546
 547         fprintf(out, "%5d%-5.5s%5.5s%5d",
 548                 resnr%100000, resname, atomname, (i+1)%100000);
 549         /* next fprintf uses built format string */
 550         if (v)
 551         {
 552             fprintf(out, format,
 553                     x[i][XX], x[i][YY], x[i][ZZ], v[i][XX], v[i][YY], v[i][ZZ]);
 554         }
 555         else
 556         {
 557             fprintf(out, format,
 558                     x[i][XX], x[i][YY], x[i][ZZ]);
 559         }
 560     }
 561
 562     write_hconf_box(out, pr, box);
 563
 564     fflush(out);
 565 }
 566
 567 void write_hconf_p(FILE *out, const char *title, t_atoms *atoms, int pr,
 568                    rvec *x, rvec *v, matrix box)
 569 {
 570     atom_id *aa;
 571     int      i;
 572
 573     snew(aa, atoms->nr);
 574     for (i = 0; (i < atoms->nr); i++)
 575     {
 576         aa[i] = i;
 577     }
 578     write_hconf_indexed_p(out, title, atoms, atoms->nr, aa, pr, x, v, box);
 579     sfree(aa);
 580 }
 581
 582 void write_conf_p(const char *outfile, const char *title,
 583                   t_atoms *atoms, int pr,
 584                   rvec *x, rvec *v, matrix box)
 585 {
 586     FILE *out;
 587
 588     out = gmx_fio_fopen(outfile, "w");
 589     write_hconf_p(out, title, atoms, pr, x, v, box);
 590     gmx_fio_fclose(out);
 591 }