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