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