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,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 "g96io.h"
  40
  41 #include <cstdio>
  42 #include <cstring>
  43
  44 #include "gromacs/fileio/trxio.h"
  45 #include "gromacs/math/vec.h"
  46 #include "gromacs/topology/atoms.h"
  47 #include "gromacs/topology/symtab.h"
  48 #include "gromacs/trajectory/trajectoryframe.h"
  49 #include "gromacs/utility/cstringutil.h"
  50 #include "gromacs/utility/fatalerror.h"
  51 #include "gromacs/utility/gmxassert.h"
  52 #include "gromacs/utility/smalloc.h"
  53
  54 #define CHAR_SHIFT 24
  55
  56 static int read_g96_pos(char line[], t_symtab *symtab,
  57                         FILE *fp, const char *infile,
  58                         t_trxframe *fr)
  59 {
  60     t_atoms   *atoms;
  61     gmx_bool   bEnd;
  62     int        nwanted, natoms, atnr, resnr = 0, oldres, newres, shift;
  63     char       anm[STRLEN], resnm[STRLEN];
  64     char       c1, c2;
  65     double     db1, db2, db3;
  66
  67     nwanted = fr->natoms;
  68
  69     if (fr->atoms != nullptr)
  70     {
  71         GMX_RELEASE_ASSERT(symtab != nullptr, "Reading a conformation from a g96 format with atom data requires a valid symbol table");
  72     }
  73     atoms = fr->atoms;
  74     if (atoms != nullptr)
  75     {
  76         atoms->haveMass    = FALSE;
  77         atoms->haveCharge  = FALSE;
  78         atoms->haveType    = FALSE;
  79         atoms->haveBState  = FALSE;
  80         atoms->havePdbInfo = FALSE;
  81     }
  82
  83     natoms = 0;
  84
  85     if (fr->bX)
  86     {
  87         if (fr->bAtoms)
  88         {
  89             shift = CHAR_SHIFT;
  90         }
  91         else
  92         {
  93             shift = 0;
  94         }
  95         newres  = -1;
  96         oldres  = -666; /* Unlikely number for the first residue! */
  97         bEnd    = FALSE;
  98         while (!bEnd && fgets2(line, STRLEN, fp))
  99         {
 100             bEnd = (std::strncmp(line, "END", 3) == 0);
 101             if (!bEnd  && (line[0] != '#'))
 102             {
 103                 if (sscanf(line+shift, "%15lf%15lf%15lf", &db1, &db2, &db3) != 3)
 104                 {
 105                     gmx_fatal(FARGS, "Did not find 3 coordinates for atom %d in %s\n",
 106                               natoms+1, infile);
 107                 }
 108                 if ((nwanted != -1) && (natoms >= nwanted))
 109                 {
 110                     gmx_fatal(FARGS,
 111                               "Found more coordinates (%d) in %s than expected %d\n",
 112                               natoms, infile, nwanted);
 113                 }
 114                 if (atoms)
 115                 {
 116                     if (fr->bAtoms &&
 117                         (sscanf(line, "%5d%c%5s%c%5s%7d", &resnr, &c1, resnm, &c2, anm, &atnr)
 118                          != 6))
 119                     {
 120                         if (oldres >= 0)
 121                         {
 122                             resnr = oldres;
 123                         }
 124                         else
 125                         {
 126                             resnr    = 1;
 127                             strncpy(resnm, "???", sizeof(resnm)-1);
 128                         }
 129                         strncpy(anm, "???", sizeof(anm)-1);
 130                     }
 131                     atoms->atomname[natoms] = put_symtab(symtab, anm);
 132                     if (resnr != oldres)
 133                     {
 134                         oldres = resnr;
 135                         newres++;
 136                         if (newres >= atoms->nr)
 137                         {
 138                             gmx_fatal(FARGS, "More residues than atoms in %s (natoms = %d)",
 139                                       infile, atoms->nr);
 140                         }
 141                         atoms->atom[natoms].resind = newres;
 142                         if (newres+1 > atoms->nres)
 143                         {
 144                             atoms->nres = newres+1;
 145                         }
 146                         t_atoms_set_resinfo(atoms, natoms, symtab, resnm, resnr, ' ', 0, ' ');
 147                     }
 148                     else
 149                     {
 150                         atoms->atom[natoms].resind = newres;
 151                     }
 152                 }
 153                 if (fr->x)
 154                 {
 155                     fr->x[natoms][0] = db1;
 156                     fr->x[natoms][1] = db2;
 157                     fr->x[natoms][2] = db3;
 158                 }
 159                 natoms++;
 160             }
 161         }
 162         if ((nwanted != -1) && natoms != nwanted)
 163         {
 164             fprintf(stderr,
 165                     "Warning: found less coordinates (%d) in %s than expected %d\n",
 166                     natoms, infile, nwanted);
 167         }
 168     }
 169
 170     fr->natoms = natoms;
 171
 172     return natoms;
 173 }
 174
 175 static int read_g96_vel(char line[], FILE *fp, const char *infile,
 176                         t_trxframe *fr)
 177 {
 178     gmx_bool   bEnd;
 179     int        nwanted, natoms = -1, shift;
 180     double     db1, db2, db3;
 181
 182     nwanted = fr->natoms;
 183
 184     if (fr->v && fr->bV)
 185     {
 186         if (strcmp(line, "VELOCITYRED") == 0)
 187         {
 188             shift = 0;
 189         }
 190         else
 191         {
 192             shift = CHAR_SHIFT;
 193         }
 194         natoms = 0;
 195         bEnd   = FALSE;
 196         while (!bEnd && fgets2(line, STRLEN, fp))
 197         {
 198             bEnd = (strncmp(line, "END", 3) == 0);
 199             if (!bEnd && (line[0] != '#'))
 200             {
 201                 if (sscanf(line+shift, "%15lf%15lf%15lf", &db1, &db2, &db3) != 3)
 202                 {
 203                     gmx_fatal(FARGS, "Did not find 3 velocities for atom %d in %s",
 204                               natoms+1, infile);
 205                 }
 206                 if ((nwanted != -1) && (natoms >= nwanted))
 207                 {
 208                     gmx_fatal(FARGS, "Found more velocities (%d) in %s than expected %d\n",
 209                               natoms, infile, nwanted);
 210                 }
 211                 if (fr->v)
 212                 {
 213                     fr->v[natoms][0] = db1;
 214                     fr->v[natoms][1] = db2;
 215                     fr->v[natoms][2] = db3;
 216                 }
 217                 natoms++;
 218             }
 219         }
 220         if ((nwanted != -1) && (natoms != nwanted))
 221         {
 222             fprintf(stderr,
 223                     "Warning: found less velocities (%d) in %s than expected %d\n",
 224                     natoms, infile, nwanted);
 225         }
 226     }
 227
 228     return natoms;
 229 }
 230
 231 int read_g96_conf(FILE *fp, const char *infile, char **name, t_trxframe *fr,
 232                   t_symtab *symtab, char *line)
 233 {
 234     gmx_bool   bAtStart, bTime, bAtoms, bPos, bVel, bBox, bEnd, bFinished;
 235     int        natoms, nbp;
 236     double     db1, db2, db3, db4, db5, db6, db7, db8, db9;
 237
 238     bAtStart = (ftell(fp) == 0);
 239
 240     clear_trxframe(fr, FALSE);
 241
 242     natoms = 0;
 243
 244     if (bAtStart)
 245     {
 246         bool foundTitle = false;
 247         while (!foundTitle && fgets2(line, STRLEN, fp))
 248         {
 249             foundTitle = (std::strcmp(line, "TITLE") == 0);
 250         }
 251         fgets2(line, STRLEN, fp);
 252         if (name != nullptr)
 253         {
 254             *name = gmx_strdup(line);
 255         }
 256         bEnd = FALSE;
 257         while (!bEnd && fgets2(line, STRLEN, fp))
 258         {
 259             bEnd = (std::strcmp(line, "END") == 0);
 260         }
 261         fgets2(line, STRLEN, fp);
 262     }
 263
 264     /* Do not get a line if we are not at the start of the file, *
 265      * because without a parameter file we don't know what is in *
 266      * the trajectory and we have already read the line in the   *
 267      * previous call (VERY DIRTY).                               */
 268     bFinished = FALSE;
 269     do
 270     {
 271         bTime  = (std::strcmp(line, "TIMESTEP") == 0);
 272         bAtoms = (std::strcmp(line, "POSITION") == 0);
 273         bPos   = (bAtoms || (strcmp(line, "POSITIONRED") == 0));
 274         bVel   = (std::strncmp(line, "VELOCITY", 8) == 0);
 275         bBox   = (std::strcmp(line, "BOX") == 0);
 276         if (bTime)
 277         {
 278             if (!fr->bTime && !fr->bX)
 279             {
 280                 fr->bStep = bTime;
 281                 fr->bTime = bTime;
 282                 do
 283                 {
 284                     bFinished = (fgets2(line, STRLEN, fp) == nullptr);
 285                 }
 286                 while (!bFinished && (line[0] == '#'));
 287                 sscanf(line, "%15" SCNd64 "%15lf", &(fr->step), &db1);
 288                 fr->time = db1;
 289             }
 290             else
 291             {
 292                 bFinished = TRUE;
 293             }
 294         }
 295         if (bPos)
 296         {
 297             if (!fr->bX)
 298             {
 299                 fr->bAtoms = bAtoms;
 300                 fr->bX     = bPos;
 301                 natoms     = read_g96_pos(line, symtab, fp, infile, fr);
 302             }
 303             else
 304             {
 305                 bFinished = TRUE;
 306             }
 307         }
 308         if (fr->v && bVel)
 309         {
 310             fr->bV = bVel;
 311             natoms = read_g96_vel(line, fp, infile, fr);
 312         }
 313         if (bBox)
 314         {
 315             fr->bBox = bBox;
 316             clear_mat(fr->box);
 317             bEnd = FALSE;
 318             while (!bEnd && fgets2(line, STRLEN, fp))
 319             {
 320                 bEnd = (strncmp(line, "END", 3) == 0);
 321                 if (!bEnd && (line[0] != '#'))
 322                 {
 323                     nbp = sscanf(line, "%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf",
 324                                  &db1, &db2, &db3, &db4, &db5, &db6, &db7, &db8, &db9);
 325                     if (nbp < 3)
 326                     {
 327                         gmx_fatal(FARGS, "Found a BOX line, but no box in %s", infile);
 328                     }
 329                     fr->box[XX][XX] = db1;
 330                     fr->box[YY][YY] = db2;
 331                     fr->box[ZZ][ZZ] = db3;
 332                     if (nbp == 9)
 333                     {
 334                         fr->box[XX][YY] = db4;
 335                         fr->box[XX][ZZ] = db5;
 336                         fr->box[YY][XX] = db6;
 337                         fr->box[YY][ZZ] = db7;
 338                         fr->box[ZZ][XX] = db8;
 339                         fr->box[ZZ][YY] = db9;
 340                     }
 341                 }
 342             }
 343             bFinished = TRUE;
 344         }
 345     }
 346     while (!bFinished && (fgets2(line, STRLEN, fp) != nullptr));
 347
 348     fr->natoms = natoms;
 349
 350     return natoms;
 351 }
 352
 353 void write_g96_conf(FILE *out, const char *title, const t_trxframe *fr,
 354                     int nindex, const int *index)
 355 {
 356     t_atoms *atoms;
 357     int      nout, i, a;
 358
 359     atoms = fr->atoms;
 360
 361     if (index)
 362     {
 363         nout = nindex;
 364     }
 365     else
 366     {
 367         nout = fr->natoms;
 368     }
 369
 370     fprintf(out, "TITLE\n%s\nEND\n", title);
 371     if (fr->bStep || fr->bTime)
 372     {
 373         /* Officially the time format is %15.9, which is not enough for 10 ns */
 374         fprintf(out, "TIMESTEP\n%15" PRId64 "%15.6f\nEND\n", fr->step, fr->time);
 375     }
 376     if (fr->bX)
 377     {
 378         if (fr->bAtoms)
 379         {
 380             fprintf(out, "POSITION\n");
 381             for (i = 0; i < nout; i++)
 382             {
 383                 if (index)
 384                 {
 385                     a = index[i];
 386                 }
 387                 else
 388                 {
 389                     a = i;
 390                 }
 391                 fprintf(out, "%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
 392                         (atoms->resinfo[atoms->atom[a].resind].nr) % 100000,
 393                         *atoms->resinfo[atoms->atom[a].resind].name,
 394                         *atoms->atomname[a], (i+1) % 10000000,
 395                         fr->x[a][XX], fr->x[a][YY], fr->x[a][ZZ]);
 396             }
 397         }
 398         else
 399         {
 400             fprintf(out, "POSITIONRED\n");
 401             for (i = 0; i < nout; i++)
 402             {
 403                 if (index)
 404                 {
 405                     a = index[i];
 406                 }
 407                 else
 408                 {
 409                     a = i;
 410                 }
 411                 fprintf(out, "%15.9f%15.9f%15.9f\n",
 412                         fr->x[a][XX], fr->x[a][YY], fr->x[a][ZZ]);
 413             }
 414         }
 415         fprintf(out, "END\n");
 416     }
 417     if (fr->bV)
 418     {
 419         if (fr->bAtoms)
 420         {
 421             fprintf(out, "VELOCITY\n");
 422             for (i = 0; i < nout; i++)
 423             {
 424                 if (index)
 425                 {
 426                     a = index[i];
 427                 }
 428                 else
 429                 {
 430                     a = i;
 431                 }
 432                 fprintf(out, "%5d %-5s %-5s%7d%15.9f%15.9f%15.9f\n",
 433                         (atoms->resinfo[atoms->atom[a].resind].nr) % 100000,
 434                         *atoms->resinfo[atoms->atom[a].resind].name,
 435                         *atoms->atomname[a], (i+1) % 10000000,
 436                         fr->v[a][XX], fr->v[a][YY], fr->v[a][ZZ]);
 437             }
 438         }
 439         else
 440         {
 441             fprintf(out, "VELOCITYRED\n");
 442             for (i = 0; i < nout; i++)
 443             {
 444                 if (index)
 445                 {
 446                     a = index[i];
 447                 }
 448                 else
 449                 {
 450                     a = i;
 451                 }
 452                 fprintf(out, "%15.9f%15.9f%15.9f\n",
 453                         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",
 462                 fr->box[XX][XX], fr->box[YY][YY], fr->box[ZZ][ZZ]);
 463         if ((fr->box[XX][YY] != 0.0F) || (fr->box[XX][ZZ] != 0.0F) || (fr->box[YY][XX] != 0.0F) ||
 464             (fr->box[YY][ZZ] != 0.0F) || (fr->box[ZZ][XX] != 0.0F) || (fr->box[ZZ][YY] != 0.0F))
 465         {
 466             fprintf(out, "%15.9f%15.9f%15.9f%15.9f%15.9f%15.9f",
 467                     fr->box[XX][YY], fr->box[XX][ZZ], fr->box[YY][XX],
 468                     fr->box[YY][ZZ], fr->box[ZZ][XX], fr->box[ZZ][YY]);
 469         }
 470         fprintf(out, "\n");
 471         fprintf(out, "END\n");
 472     }
 473 }