src/gromacs/fileio/tngio_for_tools.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include "tngio_for_tools.h"
  36
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #include "tngio.h"
  42 #include "trx.h"
  43
  44 #ifdef GMX_USE_TNG
  45 #include "../../external/tng_io/include/tng_io.h"
  46 #endif
  47
  48 #include "gromacs/utility/common.h"
  49 #include "gromacs/legacyheaders/types/atoms.h"
  50 #include "gromacs/legacyheaders/smalloc.h"
  51 #include "gromacs/legacyheaders/physics.h"
  52 #include "gromacs/legacyheaders/gmx_fatal.h"
  53
  54 void gmx_prepare_tng_writing(const char              *filename,
  55                              char                     mode,
  56                              tng_trajectory_t        *input,
  57                              tng_trajectory_t        *output,
  58                              int                      nAtoms,
  59                              const gmx_mtop_t        *mtop,
  60                              const atom_id           *index,
  61                              const char              *indexGroupName)
  62 {
  63 #ifdef GMX_USE_TNG
  64     /* FIXME after 5.0: Currently only standard block types are read */
  65     const int           defaultNumIds              = 5;
  66     static gmx_int64_t  fallbackIds[defaultNumIds] =
  67     {
  68         TNG_TRAJ_BOX_SHAPE, TNG_TRAJ_POSITIONS,
  69         TNG_TRAJ_VELOCITIES, TNG_TRAJ_FORCES,
  70         TNG_GMX_LAMBDA
  71     };
  72     static char         fallbackNames[defaultNumIds][32] =
  73     {
  74         "BOX SHAPE", "POSITIONS", "VELOCITIES",
  75         "FORCES", "LAMBDAS"
  76     };
  77
  78
  79     gmx_tng_open(filename, mode, output);
  80
  81     /* Do we have an input file in TNG format? If so, then there's
  82        more data we can copy over, rather than having to improvise. */
  83     if (*input)
  84     {
  85         /* Set parameters (compression, time per frame, molecule
  86          * information, number of frames per frame set and writing
  87          * intervals of positions, box shape and lambdas) of the
  88          * output tng container based on their respective values int
  89          * the input tng container */
  90         double      time, compression_precision;
  91         gmx_int64_t n_frames_per_frame_set, interval = -1;
  92
  93         tng_compression_precision_get(*input, &compression_precision);
  94         tng_compression_precision_set(*output, compression_precision);
  95         // TODO make this configurable in a future version
  96         char compression_type = TNG_TNG_COMPRESSION;
  97
  98         tng_molecule_system_copy(*input, *output);
  99
 100         tng_time_per_frame_get(*input, &time);
 101         tng_time_per_frame_set(*output, time);
 102
 103         tng_num_frames_per_frame_set_get(*input, &n_frames_per_frame_set);
 104         tng_num_frames_per_frame_set_set(*output, n_frames_per_frame_set);
 105
 106         for (int i = 0; i < defaultNumIds; i++)
 107         {
 108             if (tng_data_get_stride_length(*input, fallbackIds[i], -1, &interval)
 109                 == TNG_SUCCESS)
 110             {
 111                 switch (fallbackIds[i])
 112                 {
 113                     case TNG_TRAJ_POSITIONS:
 114                     case TNG_TRAJ_VELOCITIES:
 115                         tng_util_generic_write_interval_set(*output, interval, 3, fallbackIds[i],
 116                                                             fallbackNames[i], TNG_PARTICLE_BLOCK_DATA,
 117                                                             compression_type);
 118                         break;
 119                     case TNG_TRAJ_FORCES:
 120                         tng_util_generic_write_interval_set(*output, interval, 3, fallbackIds[i],
 121                                                             fallbackNames[i], TNG_PARTICLE_BLOCK_DATA,
 122                                                             TNG_GZIP_COMPRESSION);
 123                         break;
 124                     case TNG_TRAJ_BOX_SHAPE:
 125                         tng_util_generic_write_interval_set(*output, interval, 9, fallbackIds[i],
 126                                                             fallbackNames[i], TNG_NON_PARTICLE_BLOCK_DATA,
 127                                                             TNG_GZIP_COMPRESSION);
 128                         break;
 129                     case TNG_GMX_LAMBDA:
 130                         tng_util_generic_write_interval_set(*output, interval, 1, fallbackIds[i],
 131                                                             fallbackNames[i], TNG_NON_PARTICLE_BLOCK_DATA,
 132                                                             TNG_GZIP_COMPRESSION);
 133                     default:
 134                         continue;
 135                 }
 136             }
 137         }
 138
 139     }
 140     else
 141     {
 142         /* TODO after trjconv is modularized: fix this so the user can
 143            change precision when they are doing an operation where
 144            this makes sense, and not otherwise.
 145
 146            char compression = bUseLossyCompression ? TNG_TNG_COMPRESSION : TNG_GZIP_COMPRESSION;
 147            gmx_tng_set_compression_precision(*output, ndec2prec(nDecimalsOfPrecision));
 148          */
 149         gmx_tng_add_mtop(*output, mtop);
 150         tng_num_frames_per_frame_set_set(*output, 1);
 151     }
 152
 153     if (index && nAtoms > 0)
 154     {
 155         gmx_tng_setup_atom_subgroup(*output, nAtoms, index, indexGroupName);
 156     }
 157
 158     /* If for some reason there are more requested atoms than there are atoms in the
 159      * molecular system create a number of implicit atoms (without atom data) to
 160      * compensate for that. */
 161     if (nAtoms >= 0)
 162     {
 163         tng_implicit_num_particles_set(*output, nAtoms);
 164     }
 165 #else
 166     GMX_UNUSED_VALUE(filename);
 167     GMX_UNUSED_VALUE(mode);
 168     GMX_UNUSED_VALUE(input);
 169     GMX_UNUSED_VALUE(output);
 170     GMX_UNUSED_VALUE(nAtoms);
 171 #endif
 172 }
 173
 174 void gmx_write_tng_from_trxframe(tng_trajectory_t        output,
 175                                  t_trxframe             *frame,
 176                                  int                     natoms)
 177 {
 178 #ifdef GMX_USE_TNG
 179     if (frame->step > 0)
 180     {
 181         double timePerFrame = frame->time * PICO / frame->step;
 182         tng_time_per_frame_set(output, timePerFrame);
 183     }
 184     if (natoms < 0)
 185     {
 186         natoms = frame->natoms;
 187     }
 188     gmx_fwrite_tng(output,
 189                    TRUE,
 190                    frame->step,
 191                    frame->time,
 192                    0,
 193                    (const rvec *) frame->box,
 194                    natoms,
 195                    (const rvec *) frame->x,
 196                    (const rvec *) frame->v,
 197                    (const rvec *) frame->f);
 198 #else
 199     GMX_UNUSED_VALUE(output);
 200     GMX_UNUSED_VALUE(frame);
 201 #endif
 202 }
 203
 204 #ifdef GMX_USE_TNG
 205 static void
 206 convert_array_to_real_array(void       *from,
 207                             real       *to,
 208                             const float fact,
 209                             const int   nAtoms,
 210                             const int   nValues,
 211                             const char  datatype)
 212 {
 213     int i, j;
 214
 215     switch (datatype)
 216     {
 217         case TNG_FLOAT_DATA:
 218             if (sizeof(real) == sizeof(float))
 219             {
 220                 if (fact == 1)
 221                 {
 222                     memcpy(to, from, nValues * sizeof(real) * nAtoms);
 223                 }
 224                 else
 225                 {
 226                     for (i = 0; i < nAtoms; i++)
 227                     {
 228                         for (j = 0; j < nValues; j++)
 229                         {
 230                             to[i*nValues+j] = (real)((float *)from)[i*nValues+j] * fact;
 231                         }
 232                     }
 233                 }
 234             }
 235             else
 236             {
 237                 for (i = 0; i < nAtoms; i++)
 238                 {
 239                     for (j = 0; j < nValues; j++)
 240                     {
 241                         to[i*nValues+j] = (real)((float *)from)[i*nValues+j] * fact;
 242                     }
 243                 }
 244             }
 245             break;
 246         case TNG_INT_DATA:
 247             for (i = 0; i < nAtoms; i++)
 248             {
 249                 for (j = 0; j < nValues; j++)
 250                 {
 251                     to[i*nValues+j] = (real)((gmx_int64_t *)from)[i*nValues+j] * fact;
 252                 }
 253             }
 254             break;
 255         case TNG_DOUBLE_DATA:
 256             if (sizeof(real) == sizeof(double))
 257             {
 258                 if (fact == 1)
 259                 {
 260                     memcpy(to, from, nValues * sizeof(real) * nAtoms);
 261                 }
 262                 else
 263                 {
 264                     for (i = 0; i < nAtoms; i++)
 265                     {
 266                         for (j = 0; j < nValues; j++)
 267                         {
 268                             to[i*nValues+j] = (real)((double *)from)[i*nValues+j] * fact;
 269                         }
 270                     }
 271                 }
 272             }
 273             else
 274             {
 275                 for (i = 0; i < nAtoms; i++)
 276                 {
 277                     for (j = 0; j < nValues; j++)
 278                     {
 279                         to[i*nValues+j] = (real)((double *)from)[i*nValues+j] * fact;
 280                     }
 281                 }
 282             }
 283             break;
 284         default:
 285             gmx_incons("Illegal datatype when converting values to a real array!");
 286             return;
 287     }
 288     return;
 289 }
 290
 291 static real getDistanceScaleFactor(tng_trajectory_t in)
 292 {
 293     gmx_int64_t exp = -1;
 294     real        distanceScaleFactor;
 295
 296     // TODO Hopefully, TNG 2.0 will do this kind of thing for us
 297     tng_distance_unit_exponential_get(in, &exp);
 298
 299     // GROMACS expects distances in nm
 300     switch (exp)
 301     {
 302         case 9:
 303             distanceScaleFactor = NANO/NANO;
 304             break;
 305         case 10:
 306             distanceScaleFactor = NANO/ANGSTROM;
 307             break;
 308         default:
 309             distanceScaleFactor = pow(10.0, exp + 9.0);
 310     }
 311
 312     return distanceScaleFactor;
 313 }
 314 #endif
 315
 316 void gmx_tng_setup_atom_subgroup(tng_trajectory_t tng,
 317                                  const int        nind,
 318                                  const atom_id   *ind,
 319                                  const char      *name)
 320 {
 321 #ifdef GMX_USE_TNG
 322     gmx_int64_t              nAtoms, cnt, nMols;
 323     tng_molecule_t           mol, iterMol;
 324     tng_chain_t              chain;
 325     tng_residue_t            res;
 326     tng_atom_t               atom;
 327     tng_function_status      stat;
 328
 329     tng_num_particles_get(tng, &nAtoms);
 330
 331     if (nAtoms == nind)
 332     {
 333         return;
 334     }
 335
 336     stat = tng_molecule_find(tng, name, -1, &mol);
 337     if (stat == TNG_SUCCESS)
 338     {
 339         tng_molecule_num_atoms_get(tng, mol, &nAtoms);
 340         tng_molecule_cnt_get(tng, mol, &cnt);
 341         if (nAtoms == nind)
 342         {
 343             stat = TNG_SUCCESS;
 344         }
 345         else
 346         {
 347             stat = TNG_FAILURE;
 348         }
 349     }
 350     if (stat == TNG_FAILURE)
 351     {
 352         /* The indexed atoms are added to one separate molecule. */
 353         tng_molecule_alloc(tng, &mol);
 354         tng_molecule_name_set(tng, mol, name);
 355         tng_molecule_chain_add(tng, mol, "", &chain);
 356
 357         for (int i = 0; i < nind; i++)
 358         {
 359             char        temp_name[256], temp_type[256];
 360
 361             /* Try to retrieve the residue name of the atom */
 362             stat = tng_residue_name_of_particle_nr_get(tng, ind[i], temp_name, 256);
 363             if (stat != TNG_SUCCESS)
 364             {
 365                 temp_name[0] = '\0';
 366             }
 367             /* Check if the molecule of the selection already contains this residue */
 368             if (tng_chain_residue_find(tng, chain, temp_name, -1, &res)
 369                 != TNG_SUCCESS)
 370             {
 371                 tng_chain_residue_add(tng, chain, temp_name, &res);
 372             }
 373             /* Try to find the original name and type of the atom */
 374             stat = tng_atom_name_of_particle_nr_get(tng, ind[i], temp_name, 256);
 375             if (stat != TNG_SUCCESS)
 376             {
 377                 temp_name[0] = '\0';
 378             }
 379             stat = tng_atom_type_of_particle_nr_get(tng, ind[i], temp_type, 256);
 380             if (stat != TNG_SUCCESS)
 381             {
 382                 temp_type[0] = '\0';
 383             }
 384             tng_residue_atom_w_id_add(tng, res, temp_name, temp_type, ind[i], &atom);
 385         }
 386         tng_molecule_existing_add(tng, &mol);
 387     }
 388     /* Set the count of the molecule containing the selected atoms to 1 and all
 389      * other molecules to 0 */
 390     tng_molecule_cnt_set(tng, mol, 1);
 391     tng_num_molecule_types_get(tng, &nMols);
 392     for (gmx_int64_t k = 0; k < nMols; k++)
 393     {
 394         tng_molecule_of_index_get(tng, k, &iterMol);
 395         if (iterMol == mol)
 396         {
 397             continue;
 398         }
 399         tng_molecule_cnt_set(tng, iterMol, 0);
 400     }
 401 #else
 402     GMX_UNUSED_VALUE(tng);
 403     GMX_UNUSED_VALUE(nind);
 404     GMX_UNUSED_VALUE(ind);
 405     GMX_UNUSED_VALUE(name);
 406 #endif
 407 }
 408
 409 /* TODO: If/when TNG acquires the ability to copy data blocks without
 410  * uncompressing them, then this implemenation should be reconsidered.
 411  * Ideally, gmx trjconv -f a.tng -o b.tng -b 10 -e 20 would be fast
 412  * and lose no information. */
 413 gmx_bool gmx_read_next_tng_frame(tng_trajectory_t            input,
 414                                  t_trxframe                 *fr,
 415                                  gmx_int64_t                *requestedIds,
 416                                  int                         numRequestedIds)
 417 {
 418 #ifdef GMX_USE_TNG
 419     gmx_bool                bOK = TRUE;
 420     tng_function_status     stat;
 421     gmx_int64_t             numberOfAtoms = -1, frameNumber = -1;
 422     gmx_int64_t             nBlocks, blockId, *blockIds = NULL;
 423     char                    datatype      = -1, codecId;
 424     void                   *values        = NULL;
 425     double                  frameTime     = -1.0;
 426     int                     size, blockDependency;
 427     float                   prec;
 428     const int               defaultNumIds = 5;
 429     static gmx_int64_t      fallbackRequestedIds[defaultNumIds] =
 430     {
 431         TNG_TRAJ_BOX_SHAPE, TNG_TRAJ_POSITIONS,
 432         TNG_TRAJ_VELOCITIES, TNG_TRAJ_FORCES,
 433         TNG_GMX_LAMBDA
 434     };
 435
 436
 437     fr->bStep     = FALSE;
 438     fr->bTime     = FALSE;
 439     fr->bLambda   = FALSE;
 440     fr->bAtoms    = FALSE;
 441     fr->bPrec     = FALSE;
 442     fr->bX        = FALSE;
 443     fr->bV        = FALSE;
 444     fr->bF        = FALSE;
 445     fr->bBox      = FALSE;
 446
 447     /* If no specific IDs were requested read all block types that can
 448      * currently be interpreted */
 449     if (!requestedIds || numRequestedIds == 0)
 450     {
 451         numRequestedIds = defaultNumIds;
 452         requestedIds    = fallbackRequestedIds;
 453     }
 454
 455     stat = tng_num_particles_get(input, &numberOfAtoms);
 456     if (stat != TNG_SUCCESS)
 457     {
 458         gmx_file("Cannot determine number of atoms from TNG file.");
 459     }
 460     fr->natoms = numberOfAtoms;
 461
 462     if (!gmx_get_tng_data_block_types_of_next_frame(input,
 463                                                     fr->step,
 464                                                     numRequestedIds,
 465                                                     requestedIds,
 466                                                     &frameNumber,
 467                                                     &nBlocks,
 468                                                     &blockIds))
 469     {
 470         return FALSE;
 471     }
 472
 473     if (nBlocks == 0)
 474     {
 475         return FALSE;
 476     }
 477
 478     for (gmx_int64_t i = 0; i < nBlocks; i++)
 479     {
 480         blockId = blockIds[i];
 481         tng_data_block_dependency_get(input, blockId, &blockDependency);
 482         if (blockDependency & TNG_PARTICLE_DEPENDENT)
 483         {
 484             stat = tng_util_particle_data_next_frame_read(input,
 485                                                           blockId,
 486                                                           &values,
 487                                                           &datatype,
 488                                                           &frameNumber,
 489                                                           &frameTime);
 490         }
 491         else
 492         {
 493             stat = tng_util_non_particle_data_next_frame_read(input,
 494                                                               blockId,
 495                                                               &values,
 496                                                               &datatype,
 497                                                               &frameNumber,
 498                                                               &frameTime);
 499         }
 500         if (stat == TNG_CRITICAL)
 501         {
 502             gmx_file("Cannot read positions from TNG file.");
 503             return FALSE;
 504         }
 505         else if (stat == TNG_FAILURE)
 506         {
 507             continue;
 508         }
 509         switch (blockId)
 510         {
 511             case TNG_TRAJ_BOX_SHAPE:
 512                 switch (datatype)
 513                 {
 514                     case TNG_INT_DATA:
 515                         size = sizeof(gmx_int64_t);
 516                         break;
 517                     case TNG_FLOAT_DATA:
 518                         size = sizeof(float);
 519                         break;
 520                     case TNG_DOUBLE_DATA:
 521                         size = sizeof(double);
 522                         break;
 523                     default:
 524                         size = 0; /* Just to make the compiler happy. */
 525                         gmx_incons("Illegal datatype of box shape values!");
 526                 }
 527                 for (int i = 0; i < DIM; i++)
 528                 {
 529                     convert_array_to_real_array((char *)(values) + size * i * DIM,
 530                                                 (real *) fr->box[i],
 531                                                 getDistanceScaleFactor(input),
 532                                                 1,
 533                                                 DIM,
 534                                                 datatype);
 535                 }
 536                 fr->bBox = TRUE;
 537                 break;
 538             case TNG_TRAJ_POSITIONS:
 539                 srenew(fr->x, fr->natoms);
 540                 convert_array_to_real_array(values,
 541                                             (real *) fr->x,
 542                                             getDistanceScaleFactor(input),
 543                                             fr->natoms,
 544                                             DIM,
 545                                             datatype);
 546                 fr->bX = TRUE;
 547                 tng_util_frame_current_compression_get(input, blockId, &codecId, &prec);
 548                 /* This must be updated if/when more lossy compression methods are added */
 549                 if (codecId == TNG_TNG_COMPRESSION)
 550                 {
 551                     fr->prec  = prec;
 552                     fr->bPrec = TRUE;
 553                 }
 554                 break;
 555             case TNG_TRAJ_VELOCITIES:
 556                 srenew(fr->v, fr->natoms);
 557                 convert_array_to_real_array(values,
 558                                             (real *) fr->v,
 559                                             getDistanceScaleFactor(input),
 560                                             fr->natoms,
 561                                             DIM,
 562                                             datatype);
 563                 fr->bV = TRUE;
 564                 tng_util_frame_current_compression_get(input, blockId, &codecId, &prec);
 565                 /* This must be updated if/when more lossy compression methods are added */
 566                 if (codecId == TNG_TNG_COMPRESSION)
 567                 {
 568                     fr->prec  = prec;
 569                     fr->bPrec = TRUE;
 570                 }
 571                 break;
 572             case TNG_TRAJ_FORCES:
 573                 srenew(fr->f, fr->natoms);
 574                 convert_array_to_real_array(values,
 575                                             (real *) fr->f,
 576                                             getDistanceScaleFactor(input),
 577                                             fr->natoms,
 578                                             DIM,
 579                                             datatype);
 580                 fr->bF = TRUE;
 581                 break;
 582             case TNG_GMX_LAMBDA:
 583                 switch (datatype)
 584                 {
 585                     case TNG_FLOAT_DATA:
 586                         fr->lambda = (*(float *)values);
 587                         break;
 588                     case TNG_DOUBLE_DATA:
 589                         fr->lambda = (*(double *)values);
 590                         break;
 591                     default:
 592                         gmx_incons("Illegal datatype lambda value!");
 593                 }
 594                 fr->bLambda = TRUE;
 595                 break;
 596             default:
 597                 gmx_warning("Illegal block type! Currently GROMACS tools can only handle certain data types. Skipping block.");
 598         }
 599         /* values does not have to be freed before reading next frame. It will
 600          * be reallocated if it is not NULL. */
 601     }
 602
 603     fr->step  = (int) frameNumber;
 604     fr->bStep = TRUE;
 605     // Convert the time to ps
 606     fr->time  = frameTime / PICO;
 607     fr->bTime = TRUE;
 608
 609     /* values must be freed before leaving this function */
 610     sfree(values);
 611
 612     return bOK;
 613 #else
 614     GMX_UNUSED_VALUE(input);
 615     GMX_UNUSED_VALUE(fr);
 616     GMX_UNUSED_VALUE(requestedIds);
 617     return FALSE;
 618 #endif
 619 }
 620
 621 void gmx_print_tng_molecule_system(tng_trajectory_t input,
 622                                    FILE            *stream)
 623 {
 624 #ifdef GMX_USE_TNG
 625     gmx_int64_t        nMolecules, nChains, nResidues, nAtoms, *molCntList;
 626     tng_molecule_t     molecule;
 627     tng_chain_t        chain;
 628     tng_residue_t      residue;
 629     tng_atom_t         atom;
 630     char               str[256], varNAtoms;
 631
 632     tng_num_molecule_types_get(input, &nMolecules);
 633     tng_molecule_cnt_list_get(input, &molCntList);
 634     /* Can the number of particles change in the trajectory or is it constant? */
 635     tng_num_particles_variable_get(input, &varNAtoms);
 636
 637     for (gmx_int64_t i = 0; i < nMolecules; i++)
 638     {
 639         tng_molecule_of_index_get(input, i, &molecule);
 640         tng_molecule_name_get(input, molecule, str, 256);
 641         if (varNAtoms == TNG_CONSTANT_N_ATOMS)
 642         {
 643             if ((int)molCntList[i] == 0)
 644             {
 645                 continue;
 646             }
 647             fprintf(stream, "Molecule: %s, count: %d\n", str, (int)molCntList[i]);
 648         }
 649         else
 650         {
 651             fprintf(stream, "Molecule: %s\n", str);
 652         }
 653         tng_molecule_num_chains_get(input, molecule, &nChains);
 654         if (nChains > 0)
 655         {
 656             for (gmx_int64_t j = 0; j < nChains; j++)
 657             {
 658                 tng_molecule_chain_of_index_get(input, molecule, j, &chain);
 659                 tng_chain_name_get(input, chain, str, 256);
 660                 fprintf(stream, "\tChain: %s\n", str);
 661                 tng_chain_num_residues_get(input, chain, &nResidues);
 662                 for (gmx_int64_t k = 0; k < nResidues; k++)
 663                 {
 664                     tng_chain_residue_of_index_get(input, chain, k, &residue);
 665                     tng_residue_name_get(input, residue, str, 256);
 666                     fprintf(stream, "\t\tResidue: %s\n", str);
 667                     tng_residue_num_atoms_get(input, residue, &nAtoms);
 668                     for (gmx_int64_t l = 0; l < nAtoms; l++)
 669                     {
 670                         tng_residue_atom_of_index_get(input, residue, l, &atom);
 671                         tng_atom_name_get(input, atom, str, 256);
 672                         fprintf(stream, "\t\t\tAtom: %s", str);
 673                         tng_atom_type_get(input, atom, str, 256);
 674                         fprintf(stream, " (%s)\n", str);
 675                     }
 676                 }
 677             }
 678         }
 679         /* It is possible to have a molecule without chains, in which case
 680          * residues in the molecule can be iterated through without going
 681          * through chains. */
 682         else
 683         {
 684             tng_molecule_num_residues_get(input, molecule, &nResidues);
 685             if (nResidues > 0)
 686             {
 687                 for (gmx_int64_t k = 0; k < nResidues; k++)
 688                 {
 689                     tng_molecule_residue_of_index_get(input, molecule, k, &residue);
 690                     tng_residue_name_get(input, residue, str, 256);
 691                     fprintf(stream, "\t\tResidue: %s\n", str);
 692                     tng_residue_num_atoms_get(input, residue, &nAtoms);
 693                     for (gmx_int64_t l = 0; l < nAtoms; l++)
 694                     {
 695                         tng_residue_atom_of_index_get(input, residue, l, &atom);
 696                         tng_atom_name_get(input, atom, str, 256);
 697                         fprintf(stream, "\t\t\tAtom: %s", str);
 698                         tng_atom_type_get(input, atom, str, 256);
 699                         fprintf(stream, " (%s)\n", str);
 700                     }
 701                 }
 702             }
 703             else
 704             {
 705                 tng_molecule_num_atoms_get(input, molecule, &nAtoms);
 706                 for (gmx_int64_t l = 0; l < nAtoms; l++)
 707                 {
 708                     tng_molecule_atom_of_index_get(input, molecule, l, &atom);
 709                     tng_atom_name_get(input, atom, str, 256);
 710                     fprintf(stream, "\t\t\tAtom: %s", str);
 711                     tng_atom_type_get(input, atom, str, 256);
 712                     fprintf(stream, " (%s)\n", str);
 713                 }
 714             }
 715         }
 716     }
 717 #else
 718     GMX_UNUSED_VALUE(input);
 719     GMX_UNUSED_VALUE(stream);
 720 #endif
 721 }
 722
 723 gmx_bool gmx_get_tng_data_block_types_of_next_frame(tng_trajectory_t     input,
 724                                                     int                  frame,
 725                                                     int                  nRequestedIds,
 726                                                     gmx_int64_t         *requestedIds,
 727                                                     gmx_int64_t         *nextFrame,
 728                                                     gmx_int64_t         *nBlocks,
 729                                                     gmx_int64_t        **blockIds)
 730 {
 731 #ifdef GMX_USE_TNG
 732     tng_function_status stat;
 733
 734     stat = tng_util_trajectory_next_frame_present_data_blocks_find(input, frame,
 735                                                                    nRequestedIds, requestedIds,
 736                                                                    nextFrame,
 737                                                                    nBlocks, blockIds);
 738
 739     if (stat == TNG_CRITICAL)
 740     {
 741         gmx_file("Cannot read TNG file. Cannot find data blocks of next frame.");
 742     }
 743     else if (stat == TNG_FAILURE)
 744     {
 745         return FALSE;
 746     }
 747     return TRUE;
 748 #else
 749     GMX_UNUSED_VALUE(input);
 750     GMX_UNUSED_VALUE(frame);
 751     GMX_UNUSED_VALUE(nRequestedIds);
 752     GMX_UNUSED_VALUE(requestedIds);
 753     GMX_UNUSED_VALUE(nextFrame);
 754     GMX_UNUSED_VALUE(nBlocks);
 755     GMX_UNUSED_VALUE(blockIds);
 756     return FALSE;
 757 #endif
 758 }
 759
 760 gmx_bool gmx_get_tng_data_next_frame_of_block_type(tng_trajectory_t     input,
 761                                                    gmx_int64_t          blockId,
 762                                                    real               **values,
 763                                                    gmx_int64_t         *frameNumber,
 764                                                    double              *frameTime,
 765                                                    gmx_int64_t         *nValuesPerFrame,
 766                                                    gmx_int64_t         *nAtoms,
 767                                                    real                *prec,
 768                                                    char                *name,
 769                                                    int                  maxLen,
 770                                                    gmx_bool            *bOK)
 771 {
 772 #ifdef GMX_USE_TNG
 773     tng_function_status stat;
 774     char                datatype = -1, codecId;
 775     int                 blockDependency;
 776     void               *data = 0;
 777     float               localPrec;
 778
 779     stat = tng_data_block_name_get(input, blockId, name, maxLen);
 780     if (stat != TNG_SUCCESS)
 781     {
 782         gmx_file("Cannot read next frame of TNG file");
 783     }
 784     stat = tng_data_block_dependency_get(input, blockId, &blockDependency);
 785     if (stat != TNG_SUCCESS)
 786     {
 787         gmx_file("Cannot read next frame of TNG file");
 788     }
 789     if (blockDependency & TNG_PARTICLE_DEPENDENT)
 790     {
 791         tng_num_particles_get(input, nAtoms);
 792         stat = tng_util_particle_data_next_frame_read(input,
 793                                                       blockId,
 794                                                       &data,
 795                                                       &datatype,
 796                                                       frameNumber,
 797                                                       frameTime);
 798     }
 799     else
 800     {
 801         *nAtoms = 1; /* There are not actually any atoms, but it is used for
 802                         allocating memory */
 803         stat    = tng_util_non_particle_data_next_frame_read(input,
 804                                                              blockId,
 805                                                              &data,
 806                                                              &datatype,
 807                                                              frameNumber,
 808                                                              frameTime);
 809     }
 810     if (stat == TNG_CRITICAL)
 811     {
 812         gmx_file("Cannot read next frame of TNG file");
 813     }
 814     if (stat == TNG_FAILURE)
 815     {
 816         *bOK = TRUE;
 817         return FALSE;
 818     }
 819
 820     stat = tng_data_block_num_values_per_frame_get(input, blockId, nValuesPerFrame);
 821     if (stat != TNG_SUCCESS)
 822     {
 823         gmx_file("Cannot read next frame of TNG file");
 824     }
 825     snew(*values, sizeof(real) * *nValuesPerFrame * *nAtoms);
 826     convert_array_to_real_array(data,
 827                                 *values,
 828                                 getDistanceScaleFactor(input),
 829                                 *nAtoms,
 830                                 *nValuesPerFrame,
 831                                 datatype);
 832
 833     tng_util_frame_current_compression_get(input, blockId, &codecId, &localPrec);
 834
 835     /* This must be updated if/when more lossy compression methods are added */
 836     if (codecId != TNG_TNG_COMPRESSION)
 837     {
 838         *prec = -1.0;
 839     }
 840     else
 841     {
 842         *prec = localPrec;
 843     }
 844
 845     *bOK = TRUE;
 846     return TRUE;
 847 #else
 848     GMX_UNUSED_VALUE(input);
 849     GMX_UNUSED_VALUE(blockId);
 850     GMX_UNUSED_VALUE(values);
 851     GMX_UNUSED_VALUE(frameNumber);
 852     GMX_UNUSED_VALUE(frameTime);
 853     GMX_UNUSED_VALUE(nValuesPerFrame);
 854     GMX_UNUSED_VALUE(nAtoms);
 855     GMX_UNUSED_VALUE(prec);
 856     GMX_UNUSED_VALUE(name);
 857     GMX_UNUSED_VALUE(maxLen);
 858     GMX_UNUSED_VALUE(bOK);
 859     return FALSE;
 860 #endif
 861 }