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