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, 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/trx.h"
  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/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
  70     natoms = 0;
  71
  72     if (fr->bX)
  73     {
  74         if (fr->bAtoms)
  75         {
  76             shift = CHAR_SHIFT;
  77         }
  78         else
  79         {
  80             shift = 0;
  81         }
  82         newres  = -1;
  83         oldres  = -666; /* Unlikely number for the first residue! */
  84         bEnd    = FALSE;
  85         while (!bEnd && fgets2(line, STRLEN, fp))
  86         {
  87             bEnd = (std::strncmp(line, "END", 3) == 0);
  88             if (!bEnd  && (line[0] != '#'))
  89             {
  90                 if (sscanf(line+shift, "%15lf%15lf%15lf", &db1, &db2, &db3) != 3)
  91                 {
  92                     gmx_fatal(FARGS, "Did not find 3 coordinates for atom %d in %s\n",
  93                               natoms+1, infile);
  94                 }
  95                 if ((nwanted != -1) && (natoms >= nwanted))
  96                 {
  97                     gmx_fatal(FARGS,
  98                               "Found more coordinates (%d) in %s than expected %d\n",
  99                               natoms, infile, nwanted);
 100                 }
 101                 if (atoms)
 102                 {
 103                     if (fr->bAtoms &&
 104                         (sscanf(line, "%5d%c%5s%c%5s%7d", &resnr, &c1, resnm, &c2, anm, &atnr)
 105                          != 6))
 106                     {
 107                         if (oldres >= 0)
 108                         {
 109                             resnr = oldres;
 110                         }
 111                         else
 112                         {
 113                             resnr    = 1;
 114                             strncpy(resnm, "???", sizeof(resnm)-1);
 115                         }
 116                         strncpy(anm, "???", sizeof(anm)-1);
 117                     }
 118                     atoms->atomname[natoms] = put_symtab(symtab, anm);
 119                     if (resnr != oldres)
 120                     {
 121                         oldres = resnr;
 122                         newres++;
 123                         if (newres >= atoms->nr)
 124                         {
 125                             gmx_fatal(FARGS, "More residues than atoms in %s (natoms = %d)",
 126                                       infile, atoms->nr);
 127                         }
 128                         atoms->atom[natoms].resind = newres;
 129                         if (newres+1 > atoms->nres)
 130                         {
 131                             atoms->nres = newres+1;
 132                         }
 133                         t_atoms_set_resinfo(atoms, natoms, symtab, resnm, resnr, ' ', 0, ' ');
 134                     }
 135                     else
 136                     {
 137                         atoms->atom[natoms].resind = newres;
 138                     }
 139                 }
 140                 if (fr->x)
 141                 {
 142                     fr->x[natoms][0] = db1;
 143                     fr->x[natoms][1] = db2;
 144                     fr->x[natoms][2] = db3;
 145                 }
 146                 natoms++;
 147             }
 148         }
 149         if ((nwanted != -1) && natoms != nwanted)
 150         {
 151             fprintf(stderr,
 152                     "Warning: found less coordinates (%d) in %s than expected %d\n",
 153                     natoms, infile, nwanted);
 154         }
 155     }
 156
 157     fr->natoms = natoms;
 158
 159     return natoms;
 160 }
 161
 162 static int read_g96_vel(char line[], FILE *fp, const char *infile,
 163                         t_trxframe *fr)
 164 {
 165     gmx_bool   bEnd;
 166     int        nwanted, natoms = -1, shift;
 167     double     db1, db2, db3;
 168
 169     nwanted = fr->natoms;
 170
 171     if (fr->v && fr->bV)
 172     {
 173         if (strcmp(line, "VELOCITYRED") == 0)
 174         {
 175             shift = 0;
 176         }
 177         else
 178         {
 179             shift = CHAR_SHIFT;
 180         }
 181         natoms = 0;
 182         bEnd   = FALSE;
 183         while (!bEnd && fgets2(line, STRLEN, fp))
 184         {
 185             bEnd = (strncmp(line, "END", 3) == 0);
 186             if (!bEnd && (line[0] != '#'))
 187             {
 188                 if (sscanf(line+shift, "%15lf%15lf%15lf", &db1, &db2, &db3) != 3)
 189                 {
 190                     gmx_fatal(FARGS, "Did not find 3 velocities for atom %d in %s",
 191                               natoms+1, infile);
 192                 }
 193                 if ((nwanted != -1) && (natoms >= nwanted))
 194                 {
 195                     gmx_fatal(FARGS, "Found more velocities (%d) in %s than expected %d\n",
 196                               natoms, infile, nwanted);
 197                 }
 198                 if (fr->v)
 199                 {
 200                     fr->v[natoms][0] = db1;
 201                     fr->v[natoms][1] = db2;
 202                     fr->v[natoms][2] = db3;
 203                 }
 204                 natoms++;
 205             }
 206         }
 207         if ((nwanted != -1) && (natoms != nwanted))
 208         {
 209             fprintf(stderr,
 210                     "Warning: found less velocities (%d) in %s than expected %d\n",
 211                     natoms, infile, nwanted);
 212         }
 213     }
 214
 215     return natoms;
 216 }
 217
 218 int read_g96_conf(FILE *fp, const char *infile, t_trxframe *fr, char *line)
 219 {
 220     t_symtab  *symtab = NULL;
 221     gmx_bool   bAtStart, bTime, bAtoms, bPos, bVel, bBox, bEnd, bFinished;
 222     int        natoms, nbp;
 223     double     db1, db2, db3, db4, db5, db6, db7, db8, db9;
 224
 225     bAtStart = (ftell(fp) == 0);
 226
 227     clear_trxframe(fr, FALSE);
 228
 229     if (!symtab)
 230     {
 231         snew(symtab, 1);
 232         open_symtab(symtab);
 233     }
 234
 235     natoms = 0;
 236
 237     if (bAtStart)
 238     {
 239         while (!fr->bTitle && fgets2(line, STRLEN, fp))
 240         {
 241             fr->bTitle = (std::strcmp(line, "TITLE") == 0);
 242         }
 243         if (fr->title == NULL)
 244         {
 245             fgets2(line, STRLEN, fp);
 246             fr->title = gmx_strdup(line);
 247         }
 248         bEnd = FALSE;
 249         while (!bEnd && fgets2(line, STRLEN, fp))
 250         {
 251             bEnd = (std::strcmp(line, "END") == 0);
 252         }
 253         fgets2(line, STRLEN, fp);
 254     }
 255
 256     /* Do not get a line if we are not at the start of the file, *
 257      * because without a parameter file we don't know what is in *
 258      * the trajectory and we have already read the line in the   *
 259      * previous call (VERY DIRTY).                               */
 260     bFinished = FALSE;
 261     do
 262     {
 263         bTime  = (std::strcmp(line, "TIMESTEP") == 0);
 264         bAtoms = (std::strcmp(line, "POSITION") == 0);
 265         bPos   = (bAtoms || (strcmp(line, "POSITIONRED") == 0));
 266         bVel   = (std::strncmp(line, "VELOCITY", 8) == 0);
 267         bBox   = (std::strcmp(line, "BOX") == 0);
 268         if (bTime)
 269         {
 270             if (!fr->bTime && !fr->bX)
 271             {
 272                 fr->bStep = bTime;
 273                 fr->bTime = bTime;
 274                 do
 275                 {
 276                     bFinished = (fgets2(line, STRLEN, fp) == NULL);
 277                 }
 278                 while (!bFinished && (line[0] == '#'));
 279                 sscanf(line, "%15d%15lf", &(fr->step), &db1);
 280                 fr->time = db1;
 281             }
 282             else
 283             {
 284                 bFinished = TRUE;
 285             }
 286         }
 287         if (bPos)
 288         {
 289             if (!fr->bX)
 290             {
 291                 fr->bAtoms = bAtoms;
 292                 fr->bX     = bPos;
 293                 natoms     = read_g96_pos(line, symtab, fp, infile, fr);
 294             }
 295             else
 296             {
 297                 bFinished = TRUE;
 298             }
 299         }
 300         if (fr->v && bVel)
 301         {
 302             fr->bV = bVel;
 303             natoms = read_g96_vel(line, fp, infile, fr);
 304         }
 305         if (bBox)
 306         {
 307             fr->bBox = bBox;
 308             clear_mat(fr->box);
 309             bEnd = FALSE;
 310             while (!bEnd && fgets2(line, STRLEN, fp))
 311             {
 312                 bEnd = (strncmp(line, "END", 3) == 0);
 313                 if (!bEnd && (line[0] != '#'))
 314                 {
 315                     nbp = sscanf(line, "%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf%15lf",
 316                                  &db1, &db2, &db3, &db4, &db5, &db6, &db7, &db8, &db9);
 317                     if (nbp < 3)
 318                     {
 319                         gmx_fatal(FARGS, "Found a BOX line, but no box in %s", infile);
 320                     }
 321                     fr->box[XX][XX] = db1;
 322                     fr->box[YY][YY] = db2;
 323                     fr->box[ZZ][ZZ] = db3;
 324                     if (nbp == 9)
 325                     {
 326                         fr->box[XX][YY] = db4;
 327                         fr->box[XX][ZZ] = db5;
 328                         fr->box[YY][XX] = db6;
 329                         fr->box[YY][ZZ] = db7;
 330                         fr->box[ZZ][XX] = db8;
 331                         fr->box[ZZ][YY] = db9;
 332                     }
 333                 }
 334             }
 335             bFinished = TRUE;
 336         }
 337     }
 338     while (!bFinished && fgets2(line, STRLEN, fp));
 339
 340     free_symtab(symtab);
 341
 342     fr->natoms = natoms;
 343
 344     return natoms;
 345 }
 346
 347 void write_g96_conf(FILE *out, t_trxframe *fr,
 348                     int nindex, const atom_id *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%15d%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 }