src/gromacs/fileio/g96io.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 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "g96io.h"
  41
  42 #include <cstdio>
  43 #include <cstring>
  44
  45 #include "gromacs/fileio/trxio.h"
  46 #include "gromacs/math/vec.h"
  47 #include "gromacs/topology/atoms.h"
  48 #include "gromacs/topology/symtab.h"
  49 #include "gromacs/trajectory/trajectoryframe.h"
  50 #include "gromacs/utility/cstringutil.h"
  51 #include "gromacs/utility/fatalerror.h"
  52 #include "gromacs/utility/gmxassert.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 #define CHAR_SHIFT 24
  56
  57 static int read_g96_pos(char line[], t_symtab* symtab, FILE* fp, const char* infile, t_trxframe* fr)
  58 {
  59     t_atoms* atoms;
  60     gmx_bool bEnd;
  61     int      nwanted, natoms, atnr, resnr = 0, oldres, newres, shift;
  62     char     anm[STRLEN], resnm[STRLEN];
  63     char     c1, c2;
  64     double   db1, db2, db3;
  65
  66     nwanted = fr->natoms;
  67
  68     if (fr->atoms != nullptr)
  69     {
  70         GMX_RELEASE_ASSERT(symtab != nullptr,
  71                            "Reading a conformation from a g96 format with atom data requires a "
  72                            "valid symbol table");
  73     }
  74     atoms = fr->atoms;
  75     if (atoms != nullptr)
  76     {
  77         atoms->haveMass    = FALSE;
  78         atoms->haveCharge  = FALSE;
  79         atoms->haveType    = FALSE;
  80         atoms->haveBState  = FALSE;
  81         atoms->havePdbInfo = FALSE;
  82     }
  83
  84     natoms = 0;
  85
  86     if (fr->bX)
  87     {
  88         if (fr->bAtoms)
  89         {
  90             shift = CHAR_SHIFT;
  91         }
  92         else
  93         {
  94             shift = 0;
  95         }
  96         newres = -1;
  97         oldres = -666; /* Unlikely number for the first residue! */
  98         bEnd   = FALSE;
  99         while (!bEnd && fgets2(line, STRLEN, fp))
 100         {
 101             bEnd = (std::strncmp(line, "END", 3) == 0);
 102             if (!bEnd && (line[0] != '#'))
 103             {
 104                 if (sscanf(line + shift, "%15lf%15lf%15lf", &db1, &db2, &db3) != 3)
 105                 {
 106                     gmx_fatal(FARGS, "Did not find 3 coordinates for atom %d in %s\n", natoms + 1, infile);
 107                 }
 108                 if ((nwanted != -1) && (natoms >= nwanted))
 109                 {
 110                     gmx_fatal(FARGS, "Found more coordinates (%d) in %s than expected %d\n", natoms, infile, nwanted);
 111                 }
 112                 if (atoms)
 113                 {
 114                     if (fr->bAtoms
 115                         && (sscanf(line, "%5d%c%5s%c%5s%7d", &resnr, &c1, resnm, &c2, anm, &atnr) != 6))
 116                     {
 117                         if (oldres >= 0)
 118                         {
 119                             resnr = oldres;
 120                         }
 121                         else
 122                         {
 123                             resnr = 1;
 124                             strncpy(resnm, "???", sizeof(resnm) - 1);
 125                         }
 126                         strncpy(anm, "???", sizeof(anm) - 1);
 127                     }
 128                     atoms->atomname[natoms] = put_symtab(symtab, anm);
 129                     if (resnr != oldres)
 130                     {
 131                         oldres = resnr;
 132                         newres++;
 133                         if (newres >= atoms->nr)
 134                         {
 135                             gmx_fatal(FARGS, "More residues than atoms in %s (natoms = %d)", infile, atoms->nr);
 136                         }
 137                         atoms->atom[natoms].resind = newres;
 138                         if (newres + 1 > atoms->nres)
 139                         {
 140                             atoms->nres = newres + 1;
 141                         }
 142                         t_atoms_set_resinfo(atoms, natoms, symtab, resnm, resnr, ' ', 0, ' ');
 143                     }
 144                     else
 145                     {
 146                         atoms->atom[natoms].resind = newres;
 147                     }
 148                 }
 149                 if (fr->x)
 150                 {
 151                     fr->x[natoms][0] = db1;
 152                     fr->x[natoms][1] = db2;
 153                     fr->x[natoms][2] = db3;
 154                 }
 155                 natoms++;
 156             }
 157         }
 158         if ((nwanted != -1) && natoms != nwanted)
 159         {
 160             fprintf(stderr, "Warning: found less coordinates (%d) in %s than expected %d\n", natoms, infile, nwanted);
 161         }
 162     }
 163
 164     fr->natoms = natoms;
 165
 166     return natoms;
 167 }
 168
 169 static int read_g96_vel(char line[], FILE* fp, const char* infile, t_trxframe* fr)
 170 {
 171     gmx_bool bEnd;
 172     int      nwanted, natoms = -1, shift;
 173     double   db1, db2, db3;
 174
 175     nwanted = fr->natoms;
 176
 177     if (fr->v && fr->bV)
 178     {
 179         if (strcmp(line, "VELOCITYRED") == 0)
 180         {
 181             shift = 0;
 182         }
 183         else
 184         {
 185             shift = CHAR_SHIFT;
 186         }
 187         natoms = 0;
 188         bEnd   = FALSE;
 189         while (!bEnd && fgets2(line, STRLEN, fp))
 190         {
 191             bEnd = (strncmp(line, "END", 3) == 0);
 192             if (!bEnd && (line[0] != '#'))
 193             {
 194                 if (sscanf(line + shift, "%15lf%15lf%15lf", &db1, &db2, &db3) != 3)
 195                 {
 196                     gmx_fatal(FARGS, "Did not find 3 velocities for atom %d in %s", natoms + 1, infile);
 197                 }
 198                 if ((nwanted != -1) && (natoms >= nwanted))
 199                 {
 200                     gmx_fatal(FARGS, "Found more velocities (%d) in %s than expected %d\n", natoms, infile, nwanted);
 201                 }
 202                 if (fr->v)
 203                 {
 204                     fr->v[natoms][0] = db1;
 205                     fr->v[natoms][1] = db2;
 206                     fr->v[natoms][2] = db3;
 207                 }
 208                 natoms++;
 209             }
 210         }
 211         if ((nwanted != -1) && (natoms != nwanted))
 212         {
 213             fprintf(stderr, "Warning: found less velocities (%d) in %s than expected %d\n", natoms, infile, nwanted);
 214         }
 215     }
 216
 217     return natoms;
 218 }
 219
 220 int read_g96_conf(FILE* fp, const char* infile, char** name, t_trxframe* fr, t_symtab* symtab, char* line)
 221 {
 222     gmx_bool bAtStart, bTime, bAtoms, bPos, bVel, bBox, bEnd, bFinished;
 223     int      natoms, nbp;
 224     double   db1, db2, db3, db4, db5, db6, db7, db8, db9;
 225
 226     bAtStart = (ftell(fp) == 0);
 227
 228     clear_trxframe(fr, FALSE);
 229
 230     natoms = 0;
 231
 232     if (bAtStart)
 233     {
 234         bool foundTitle = false;
 235         while (!foundTitle && fgets2(line, STRLEN, fp))
 236         {
 237             foundTitle = (std::strcmp(line, "TITLE") == 0);
 238         }
 239         fgets2(line, STRLEN, fp);
 240         if (name != nullptr)
 241         {
 242             *name = gmx_strdup(line);
 243         }
 244         bEnd = FALSE;
 245         while (!bEnd && fgets2(line, STRLEN, fp))
 246         {
 247             bEnd = (std::strcmp(line, "END") == 0);
 248         }
 249         fgets2(line, STRLEN, fp);
 250     }
 251
 252     /* Do not get a line if we are not at the start of the file, *
 253      * because without a parameter file we don't know what is in *
 254      * the trajectory and we have already read the line in the   *
 255      * previous call (VERY DIRTY).                               */
 256     bFinished = FALSE;
 257     do
 258     {
 259         bTime  = (std::strcmp(line, "TIMESTEP") == 0);
 260         bAtoms = (std::strcmp(line, "POSITION") == 0);
 261         bPos   = (bAtoms || (strcmp(line, "POSITIONRED") == 0));
 262         bVel   = (std::strncmp(line, "VELOCITY", 8) == 0);
 263         bBox   = (std::strcmp(line, "BOX") == 0);
 264         if (bTime)
 265         {
 266             if (!fr->bTime && !fr->bX)
 267             {
 268                 fr->bStep = bTime;
 269                 fr->bTime = bTime;
 270                 do
 271                 {
 272                     bFinished = (fgets2(line, STRLEN, fp) == nullptr);
 273                 } while (!bFinished && (line[0] == '#'));
 274                 sscanf(line, "%15" SCNd64 "%15lf", &(fr->step), &db1);
 275                 fr->time = db1;
 276             }
 277             else
 278             {
 279                 bFinished = TRUE;
 280             }
 281         }
 282         if (bPos)
 283         {
 284             if (!fr->bX)
 285             {
 286                 fr->bAtoms = bAtoms;
 287                 fr->bX     = bPos;
 288                 natoms     = read_g96_pos(line, symtab, fp, infile, fr);
 289             }
 290             else
 291             {
 292                 bFinished = TRUE;
 293             }
 294         }
 295         if (fr->v && bVel)
 296         {
 297             fr->bV = bVel;
 298             natoms = read_g96_vel(line, fp, infile, fr);
 299         }
 300         if (bBox)
 301         {
 302             fr->bBox = bBox;
 303             clear_mat(fr->box);
 304             bEnd = FALSE;
 305             while (!bEnd && fgets2(line, STRLEN, fp))
 306             {
 307                 bEnd = (strncmp(line, "END", 3) == 0);
 308                 if (!bEnd && (line[0] != '#'))
 309                 {
 310                     nbp = sscanf(line,
 311                                  "%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf",
 312                                  &db1,
 313                                  &db2,
 314                                  &db3,
 315                                  &db4,
 316                                  &db5,
 317                                  &db6,
 318                                  &db7,
 319                                  &db8,
 320                                  &db9);
 321                     if (nbp < 3)
 322                     {
 323                         gmx_fatal(FARGS, "Found a BOX line, but no box in %s", infile);
 324                     }
 325                     fr->box[XX][XX] = db1;
 326                     fr->box[YY][YY] = db2;
 327                     fr->box[ZZ][ZZ] = db3;
 328                     if (nbp == 9)
 329                     {
 330                         fr->box[XX][YY] = db4;
 331                         fr->box[XX][ZZ] = db5;
 332                         fr->box[YY][XX] = db6;
 333                         fr->box[YY][ZZ] = db7;
 334                         fr->box[ZZ][XX] = db8;
 335                         fr->box[ZZ][YY] = db9;
 336                     }
 337                 }
 338             }
 339             bFinished = TRUE;
 340         }
 341     } while (!bFinished && (fgets2(line, STRLEN, fp) != nullptr));
 342
 343     fr->natoms = natoms;
 344
 345     return natoms;
 346 }
 347
 348 void write_g96_conf(FILE* out, const char* title, const t_trxframe* fr, int nindex, const int* index)
 349 {
 350     t_atoms* atoms;
 351     int      nout, i, a;
 352
 353     atoms = fr->atoms;
 354
 355     if (index)
 356     {
 357         nout = nindex;
 358     }
 359     else
 360     {
 361         nout = fr->natoms;
 362     }
 363
 364     fprintf(out, "TITLE\n%s\nEND\n", title);
 365     if (fr->bStep || fr->bTime)
 366     {
 367         /* Officially the time format is %15.9, which is not enough for 10 ns */
 368         fprintf(out, "TIMESTEP\n%15" PRId64 "%15.6f\nEND\n", fr->step, fr->time);
 369     }
 370     if (fr->bX)
 371     {
 372         if (fr->bAtoms)
 373         {
 374             fprintf(out, "POSITION\n");
 375             for (i = 0; i < nout; i++)
 376             {
 377                 if (index)
 378                 {
 379                     a = index[i];
 380                 }
 381                 else
 382                 {
 383                     a = i;
 384                 }
 385                 fprintf(out,
 386                         "%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
 387                         (atoms->resinfo[atoms->atom[a].resind].nr) % 100000,
 388                         *atoms->resinfo[atoms->atom[a].resind].name,
 389                         *atoms->atomname[a],
 390                         (i + 1) % 10000000,
 391                         fr->x[a][XX],
 392                         fr->x[a][YY],
 393                         fr->x[a][ZZ]);
 394             }
 395         }
 396         else
 397         {
 398             fprintf(out, "POSITIONRED\n");
 399             for (i = 0; i < nout; i++)
 400             {
 401                 if (index)
 402                 {
 403                     a = index[i];
 404                 }
 405                 else
 406                 {
 407                     a = i;
 408                 }
 409                 fprintf(out, "%15.9f%15.9f%15.9f\n", fr->x[a][XX], fr->x[a][YY], fr->x[a][ZZ]);
 410             }
 411         }
 412         fprintf(out, "END\n");
 413     }
 414     if (fr->bV)
 415     {
 416         if (fr->bAtoms)
 417         {
 418             fprintf(out, "VELOCITY\n");
 419             for (i = 0; i < nout; i++)
 420             {
 421                 if (index)
 422                 {
 423                     a = index[i];
 424                 }
 425                 else
 426                 {
 427                     a = i;
 428                 }
 429                 fprintf(out,
 430                         "%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
 431                         (atoms->resinfo[atoms->atom[a].resind].nr) % 100000,
 432                         *atoms->resinfo[atoms->atom[a].resind].name,
 433                         *atoms->atomname[a],
 434                         (i + 1) % 10000000,
 435                         fr->v[a][XX],
 436                         fr->v[a][YY],
 437                         fr->v[a][ZZ]);
 438             }
 439         }
 440         else
 441         {
 442             fprintf(out, "VELOCITYRED\n");
 443             for (i = 0; i < nout; i++)
 444             {
 445                 if (index)
 446                 {
 447                     a = index[i];
 448                 }
 449                 else
 450                 {
 451                     a = i;
 452                 }
 453                 fprintf(out, "%15.9f%15.9f%15.9f\n", fr->v[a][XX], fr->v[a][YY], fr->v[a][ZZ]);
 454             }
 455         }
 456         fprintf(out, "END\n");
 457     }
 458     if (fr->bBox)
 459     {
 460         fprintf(out, "BOX\n");
 461         fprintf(out, "%15.9f%15.9f%15.9f", fr->box[XX][XX], fr->box[YY][YY], fr->box[ZZ][ZZ]);
 462         if ((fr->box[XX][YY] != 0.0F) || (fr->box[XX][ZZ] != 0.0F) || (fr->box[YY][XX] != 0.0F)
 463             || (fr->box[YY][ZZ] != 0.0F) || (fr->box[ZZ][XX] != 0.0F) || (fr->box[ZZ][YY] != 0.0F))
 464         {
 465             fprintf(out,
 466                     "%15.9f%15.9f%15.9f%15.9f%15.9f%15.9f",
 467                     fr->box[XX][YY],
 468                     fr->box[XX][ZZ],
 469                     fr->box[YY][XX],
 470                     fr->box[YY][ZZ],
 471                     fr->box[ZZ][XX],
 472                     fr->box[ZZ][YY]);
 473         }
 474         fprintf(out, "\n");
 475         fprintf(out, "END\n");
 476     }
 477 }