src/gromacs/fileio/tngio_for_tools.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014, 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 <math.h>
  42
  43 #include "tngio.h"
  44 #include "trx.h"
  45
  46 #ifdef GMX_USE_TNG
  47 #include "../../external/tng_io/include/tng_io.h"
  48 #endif
  49
  50 #include "gromacs/math/units.h"
  51 #include "gromacs/utility/common.h"
  52 #include "gromacs/utility/fatalerror.h"
  53 #include "gromacs/utility/smalloc.h"
  54
  55 void gmx_prepare_tng_writing(const char              *filename,
  56                              char                     mode,
  57                              tng_trajectory_t        *input,
  58                              tng_trajectory_t        *output,
  59                              int                      nAtoms,
  60                              const gmx_mtop_t        *mtop,
  61                              const atom_id           *index,
  62                              const char              *indexGroupName)
  63 {
  64 #ifdef GMX_USE_TNG
  65     /* FIXME after 5.0: Currently only standard block types are read */
  66     const int           defaultNumIds              = 5;
  67     static gmx_int64_t  fallbackIds[defaultNumIds] =
  68     {
  69         TNG_TRAJ_BOX_SHAPE, TNG_TRAJ_POSITIONS,
  70         TNG_TRAJ_VELOCITIES, TNG_TRAJ_FORCES,
  71         TNG_GMX_LAMBDA
  72     };
  73     static char         fallbackNames[defaultNumIds][32] =
  74     {
  75         "BOX SHAPE", "POSITIONS", "VELOCITIES",
  76         "FORCES", "LAMBDAS"
  77     };
  78
  79
  80     gmx_tng_open(filename, mode, output);
  81
  82     /* Do we have an input file in TNG format? If so, then there's
  83        more data we can copy over, rather than having to improvise. */
  84     if (*input)
  85     {
  86         /* Set parameters (compression, time per frame, molecule
  87          * information, number of frames per frame set and writing
  88          * intervals of positions, box shape and lambdas) of the
  89          * output tng container based on their respective values int
  90          * the input tng container */
  91         double      time, compression_precision;
  92         gmx_int64_t n_frames_per_frame_set, interval = -1;
  93
  94         tng_compression_precision_get(*input, &compression_precision);
  95         tng_compression_precision_set(*output, compression_precision);
  96         // TODO make this configurable in a future version
  97         char compression_type = TNG_TNG_COMPRESSION;
  98
  99         tng_molecule_system_copy(*input, *output);
 100
 101         tng_time_per_frame_get(*input, &time);
 102         tng_time_per_frame_set(*output, time);
 103
 104         tng_num_frames_per_frame_set_get(*input, &n_frames_per_frame_set);
 105         tng_num_frames_per_frame_set_set(*output, n_frames_per_frame_set);
 106
 107         for (int i = 0; i < defaultNumIds; i++)
 108         {
 109             if (tng_data_get_stride_length(*input, fallbackIds[i], -1, &interval)
 110                 == TNG_SUCCESS)
 111             {
 112                 switch (fallbackIds[i])
 113                 {
 114                     case TNG_TRAJ_POSITIONS:
 115                     case TNG_TRAJ_VELOCITIES:
 116                         tng_util_generic_write_interval_set(*output, interval, 3, fallbackIds[i],
 117                                                             fallbackNames[i], TNG_PARTICLE_BLOCK_DATA,
 118                                                             compression_type);
 119                         break;
 120                     case TNG_TRAJ_FORCES:
 121                         tng_util_generic_write_interval_set(*output, interval, 3, fallbackIds[i],
 122                                                             fallbackNames[i], TNG_PARTICLE_BLOCK_DATA,
 123                                                             TNG_GZIP_COMPRESSION);
 124                         break;
 125                     case TNG_TRAJ_BOX_SHAPE:
 126                         tng_util_generic_write_interval_set(*output, interval, 9, fallbackIds[i],
 127                                                             fallbackNames[i], TNG_NON_PARTICLE_BLOCK_DATA,
 128                                                             TNG_GZIP_COMPRESSION);
 129                         break;
 130                     case TNG_GMX_LAMBDA:
 131                         tng_util_generic_write_interval_set(*output, interval, 1, fallbackIds[i],
 132                                                             fallbackNames[i], TNG_NON_PARTICLE_BLOCK_DATA,
 133                                                             TNG_GZIP_COMPRESSION);
 134                     default:
 135                         continue;
 136                 }
 137             }
 138         }
 139
 140     }
 141     else
 142     {
 143         /* TODO after trjconv is modularized: fix this so the user can
 144            change precision when they are doing an operation where
 145            this makes sense, and not otherwise.
 146
 147            char compression = bUseLossyCompression ? TNG_TNG_COMPRESSION : TNG_GZIP_COMPRESSION;
 148            gmx_tng_set_compression_precision(*output, ndec2prec(nDecimalsOfPrecision));
 149          */
 150         gmx_tng_add_mtop(*output, mtop);
 151         tng_num_frames_per_frame_set_set(*output, 1);
 152     }
 153
 154     if (index && nAtoms > 0)
 155     {
 156         gmx_tng_setup_atom_subgroup(*output, nAtoms, index, indexGroupName);
 157     }
 158
 159     /* If for some reason there are more requested atoms than there are atoms in the
 160      * molecular system create a number of implicit atoms (without atom data) to
 161      * compensate for that. */
 162     if (nAtoms >= 0)
 163     {
 164         tng_implicit_num_particles_set(*output, nAtoms);
 165     }
 166 #else
 167     GMX_UNUSED_VALUE(filename);
 168     GMX_UNUSED_VALUE(mode);
 169     GMX_UNUSED_VALUE(input);
 170     GMX_UNUSED_VALUE(output);
 171     GMX_UNUSED_VALUE(nAtoms);
 172 #endif
 173 }
 174
 175 void gmx_write_tng_from_trxframe(tng_trajectory_t        output,
 176                                  t_trxframe             *frame,
 177                                  int                     natoms)
 178 {
 179 #ifdef GMX_USE_TNG
 180     if (frame->step > 0)
 181     {
 182         double timePerFrame = frame->time * PICO / frame->step;
 183         tng_time_per_frame_set(output, timePerFrame);
 184     }
 185     if (natoms < 0)
 186     {
 187         natoms = frame->natoms;
 188     }
 189     gmx_fwrite_tng(output,
 190                    TRUE,
 191                    frame->step,
 192                    frame->time,
 193                    0,
 194                    (const rvec *) frame->box,
 195                    natoms,
 196                    (const rvec *) frame->x,
 197                    (const rvec *) frame->v,
 198                    (const rvec *) frame->f);
 199 #else
 200     GMX_UNUSED_VALUE(output);
 201     GMX_UNUSED_VALUE(frame);
 202 #endif
 203 }
 204
 205 #ifdef GMX_USE_TNG
 206 static void
 207 convert_array_to_real_array(void       *from,
 208                             real       *to,
 209                             const float fact,
 210                             const int   nAtoms,
 211                             const int   nValues,
 212                             const char  datatype)
 213 {
 214     int i, j;
 215
 216     switch (datatype)
 217     {
 218         case TNG_FLOAT_DATA:
 219             if (sizeof(real) == sizeof(float))
 220             {
 221                 if (fact == 1)
 222                 {
 223                     memcpy(to, from, nValues * sizeof(real) * nAtoms);
 224                 }
 225                 else
 226                 {
 227                     for (i = 0; i < nAtoms; i++)
 228                     {
 229                         for (j = 0; j < nValues; j++)
 230                         {
 231                             to[i*nValues+j] = (real)((float *)from)[i*nValues+j] * fact;
 232                         }
 233                     }
 234                 }
 235             }
 236             else
 237             {
 238                 for (i = 0; i < nAtoms; i++)
 239                 {
 240                     for (j = 0; j < nValues; j++)
 241                     {
 242                         to[i*nValues+j] = (real)((float *)from)[i*nValues+j] * fact;
 243                     }
 244                 }
 245             }
 246             break;
 247         case TNG_INT_DATA:
 248             for (i = 0; i < nAtoms; i++)
 249             {
 250                 for (j = 0; j < nValues; j++)
 251                 {
 252                     to[i*nValues+j] = (real)((gmx_int64_t *)from)[i*nValues+j] * fact;
 253                 }
 254             }
 255             break;
 256         case TNG_DOUBLE_DATA:
 257             if (sizeof(real) == sizeof(double))
 258             {
 259                 if (fact == 1)
 260                 {
 261                     memcpy(to, from, nValues * sizeof(real) * nAtoms);
 262                 }
 263                 else
 264                 {
 265                     for (i = 0; i < nAtoms; i++)
 266                     {
 267                         for (j = 0; j < nValues; j++)
 268                         {
 269                             to[i*nValues+j] = (real)((double *)from)[i*nValues+j] * fact;
 270                         }
 271                     }
 272                 }
 273             }
 274             else
 275             {
 276                 for (i = 0; i < nAtoms; i++)
 277                 {
 278                     for (j = 0; j < nValues; j++)
 279                     {
 280                         to[i*nValues+j] = (real)((double *)from)[i*nValues+j] * fact;
 281                     }
 282                 }
 283             }
 284             break;
 285         default:
 286             gmx_incons("Illegal datatype when converting values to a real array!");
 287             return;
 288     }
 289     return;
 290 }
 291
 292 static real getDistanceScaleFactor(tng_trajectory_t in)
 293 {
 294     gmx_int64_t exp = -1;
 295     real        distanceScaleFactor;
 296
 297     // TODO Hopefully, TNG 2.0 will do this kind of thing for us
 298     tng_distance_unit_exponential_get(in, &exp);
 299
 300     // GROMACS expects distances in nm
 301     switch (exp)
 302     {
 303         case 9:
 304             distanceScaleFactor = NANO/NANO;
 305             break;
 306         case 10:
 307             distanceScaleFactor = NANO/ANGSTROM;
 308             break;
 309         default:
 310             distanceScaleFactor = pow(10.0, exp + 9.0);
 311     }
 312
 313     return distanceScaleFactor;
 314 }
 315 #endif
 316
 317 void gmx_tng_setup_atom_subgroup(tng_trajectory_t tng,
 318                                  const int        nind,
 319                                  const atom_id   *ind,
 320                                  const char      *name)
 321 {
 322 #ifdef GMX_USE_TNG
 323     gmx_int64_t              nAtoms, cnt, nMols;
 324     tng_molecule_t           mol, iterMol;
 325     tng_chain_t              chain;
 326     tng_residue_t            res;
 327     tng_atom_t               atom;
 328     tng_function_status      stat;
 329
 330     tng_num_particles_get(tng, &nAtoms);
 331
 332     if (nAtoms == nind)
 333     {
 334         return;
 335     }
 336
 337     stat = tng_molecule_find(tng, name, -1, &mol);
 338     if (stat == TNG_SUCCESS)
 339     {
 340         tng_molecule_num_atoms_get(tng, mol, &nAtoms);
 341         tng_molecule_cnt_get(tng, mol, &cnt);
 342         if (nAtoms == nind)
 343         {
 344             stat = TNG_SUCCESS;
 345         }
 346         else
 347         {
 348             stat = TNG_FAILURE;
 349         }
 350     }
 351     if (stat == TNG_FAILURE)
 352     {
 353         /* The indexed atoms are added to one separate molecule. */
 354         tng_molecule_alloc(tng, &mol);
 355         tng_molecule_name_set(tng, mol, name);
 356         tng_molecule_chain_add(tng, mol, "", &chain);
 357
 358         for (int i = 0; i < nind; i++)
 359         {
 360             char        temp_name[256], temp_type[256];
 361
 362             /* Try to retrieve the residue name of the atom */
 363             stat = tng_residue_name_of_particle_nr_get(tng, ind[i], temp_name, 256);
 364             if (stat != TNG_SUCCESS)
 365             {
 366                 temp_name[0] = '\0';
 367             }
 368             /* Check if the molecule of the selection already contains this residue */
 369             if (tng_chain_residue_find(tng, chain, temp_name, -1, &res)
 370                 != TNG_SUCCESS)
 371             {
 372                 tng_chain_residue_add(tng, chain, temp_name, &res);
 373             }
 374             /* Try to find the original name and type of the atom */
 375             stat = tng_atom_name_of_particle_nr_get(tng, ind[i], temp_name, 256);
 376             if (stat != TNG_SUCCESS)
 377             {
 378                 temp_name[0] = '\0';
 379             }
 380             stat = tng_atom_type_of_particle_nr_get(tng, ind[i], temp_type, 256);
 381             if (stat != TNG_SUCCESS)
 382             {
 383                 temp_type[0] = '\0';
 384             }
 385             tng_residue_atom_w_id_add(tng, res, temp_name, temp_type, ind[i], &atom);
 386         }
 387         tng_molecule_existing_add(tng, &mol);
 388     }
 389     /* Set the count of the molecule containing the selected atoms to 1 and all
 390      * other molecules to 0 */
 391     tng_molecule_cnt_set(tng, mol, 1);
 392     tng_num_molecule_types_get(tng, &nMols);
 393     for (gmx_int64_t k = 0; k < nMols; k++)
 394     {
 395         tng_molecule_of_index_get(tng, k, &iterMol);
 396         if (iterMol == mol)
 397         {
 398             continue;
 399         }
 400         tng_molecule_cnt_set(tng, iterMol, 0);
 401     }
 402 #else
 403     GMX_UNUSED_VALUE(tng);
 404     GMX_UNUSED_VALUE(nind);
 405     GMX_UNUSED_VALUE(ind);
 406     GMX_UNUSED_VALUE(name);
 407 #endif
 408 }
 409
 410 /* TODO: If/when TNG acquires the ability to copy data blocks without
 411  * uncompressing them, then this implemenation should be reconsidered.
 412  * Ideally, gmx trjconv -f a.tng -o b.tng -b 10 -e 20 would be fast
 413  * and lose no information. */
 414 gmx_bool gmx_read_next_tng_frame(tng_trajectory_t            input,
 415                                  t_trxframe                 *fr,
 416                                  gmx_int64_t                *requestedIds,
 417                                  int                         numRequestedIds)
 418 {
 419 #ifdef GMX_USE_TNG
 420     gmx_bool                bOK = TRUE;
 421     tng_function_status     stat;
 422     gmx_int64_t             numberOfAtoms = -1, frameNumber = -1;
 423     gmx_int64_t             nBlocks, blockId, *blockIds = NULL, codecId;
 424     char                    datatype      = -1;
 425     void                   *values        = NULL;
 426     double                  frameTime     = -1.0;
 427     int                     size, blockDependency;
 428     float                   prec;
 429     const int               defaultNumIds = 5;
 430     static gmx_int64_t      fallbackRequestedIds[defaultNumIds] =
 431     {
 432         TNG_TRAJ_BOX_SHAPE, TNG_TRAJ_POSITIONS,
 433         TNG_TRAJ_VELOCITIES, TNG_TRAJ_FORCES,
 434         TNG_GMX_LAMBDA
 435     };
 436
 437
 438     fr->bStep     = FALSE;
 439     fr->bTime     = FALSE;
 440     fr->bLambda   = FALSE;
 441     fr->bAtoms    = FALSE;
 442     fr->bPrec     = FALSE;
 443     fr->bX        = FALSE;
 444     fr->bV        = FALSE;
 445     fr->bF        = FALSE;
 446     fr->bBox      = FALSE;
 447
 448     /* If no specific IDs were requested read all block types that can
 449      * currently be interpreted */
 450     if (!requestedIds || numRequestedIds == 0)
 451     {
 452         numRequestedIds = defaultNumIds;
 453         requestedIds    = fallbackRequestedIds;
 454     }
 455
 456     stat = tng_num_particles_get(input, &numberOfAtoms);
 457     if (stat != TNG_SUCCESS)
 458     {
 459         gmx_file("Cannot determine number of atoms from TNG file.");
 460     }
 461     fr->natoms = numberOfAtoms;
 462
 463     if (!gmx_get_tng_data_block_types_of_next_frame(input,
 464                                                     fr->step,
 465                                                     numRequestedIds,
 466                                                     requestedIds,
 467                                                     &frameNumber,
 468                                                     &nBlocks,
 469                                                     &blockIds))
 470     {
 471         return FALSE;
 472     }
 473
 474     if (nBlocks == 0)
 475     {
 476         return FALSE;
 477     }
 478
 479     for (gmx_int64_t i = 0; i < nBlocks; i++)
 480     {
 481         blockId = blockIds[i];
 482         tng_data_block_dependency_get(input, blockId, &blockDependency);
 483         if (blockDependency & TNG_PARTICLE_DEPENDENT)
 484         {
 485             stat = tng_util_particle_data_next_frame_read(input,
 486                                                           blockId,
 487                                                           &values,
 488                                                           &datatype,
 489                                                           &frameNumber,
 490                                                           &frameTime);
 491         }
 492         else
 493         {
 494             stat = tng_util_non_particle_data_next_frame_read(input,
 495                                                               blockId,
 496                                                               &values,
 497                                                               &datatype,
 498                                                               &frameNumber,
 499                                                               &frameTime);
 500         }
 501         if (stat == TNG_CRITICAL)
 502         {
 503             gmx_file("Cannot read positions from TNG file.");
 504             return FALSE;
 505         }
 506         else if (stat == TNG_FAILURE)
 507         {
 508             continue;
 509         }
 510         switch (blockId)
 511         {
 512             case TNG_TRAJ_BOX_SHAPE:
 513                 switch (datatype)
 514                 {
 515                     case TNG_INT_DATA:
 516                         size = sizeof(gmx_int64_t);
 517                         break;
 518                     case TNG_FLOAT_DATA:
 519                         size = sizeof(float);
 520                         break;
 521                     case TNG_DOUBLE_DATA:
 522                         size = sizeof(double);
 523                         break;
 524                     default:
 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;
 775     gmx_int64_t         codecId;
 776     int                 blockDependency;
 777     void               *data = 0;
 778     float               localPrec;
 779
 780     stat = tng_data_block_name_get(input, blockId, name, maxLen);
 781     if (stat != TNG_SUCCESS)
 782     {
 783         gmx_file("Cannot read next frame of TNG file");
 784     }
 785     stat = tng_data_block_dependency_get(input, blockId, &blockDependency);
 786     if (stat != TNG_SUCCESS)
 787     {
 788         gmx_file("Cannot read next frame of TNG file");
 789     }
 790     if (blockDependency & TNG_PARTICLE_DEPENDENT)
 791     {
 792         tng_num_particles_get(input, nAtoms);
 793         stat = tng_util_particle_data_next_frame_read(input,
 794                                                       blockId,
 795                                                       &data,
 796                                                       &datatype,
 797                                                       frameNumber,
 798                                                       frameTime);
 799     }
 800     else
 801     {
 802         *nAtoms = 1; /* There are not actually any atoms, but it is used for
 803                         allocating memory */
 804         stat    = tng_util_non_particle_data_next_frame_read(input,
 805                                                              blockId,
 806                                                              &data,
 807                                                              &datatype,
 808                                                              frameNumber,
 809                                                              frameTime);
 810     }
 811     if (stat == TNG_CRITICAL)
 812     {
 813         gmx_file("Cannot read next frame of TNG file");
 814     }
 815     if (stat == TNG_FAILURE)
 816     {
 817         *bOK = TRUE;
 818         return FALSE;
 819     }
 820
 821     stat = tng_data_block_num_values_per_frame_get(input, blockId, nValuesPerFrame);
 822     if (stat != TNG_SUCCESS)
 823     {
 824         gmx_file("Cannot read next frame of TNG file");
 825     }
 826     snew(*values, sizeof(real) * *nValuesPerFrame * *nAtoms);
 827     convert_array_to_real_array(data,
 828                                 *values,
 829                                 getDistanceScaleFactor(input),
 830                                 *nAtoms,
 831                                 *nValuesPerFrame,
 832                                 datatype);
 833
 834     tng_util_frame_current_compression_get(input, blockId, &codecId, &localPrec);
 835
 836     /* This must be updated if/when more lossy compression methods are added */
 837     if (codecId != TNG_TNG_COMPRESSION)
 838     {
 839         *prec = -1.0;
 840     }
 841     else
 842     {
 843         *prec = localPrec;
 844     }
 845
 846     *bOK = TRUE;
 847     return TRUE;
 848 #else
 849     GMX_UNUSED_VALUE(input);
 850     GMX_UNUSED_VALUE(blockId);
 851     GMX_UNUSED_VALUE(values);
 852     GMX_UNUSED_VALUE(frameNumber);
 853     GMX_UNUSED_VALUE(frameTime);
 854     GMX_UNUSED_VALUE(nValuesPerFrame);
 855     GMX_UNUSED_VALUE(nAtoms);
 856     GMX_UNUSED_VALUE(prec);
 857     GMX_UNUSED_VALUE(name);
 858     GMX_UNUSED_VALUE(maxLen);
 859     GMX_UNUSED_VALUE(bOK);
 860     return FALSE;
 861 #endif
 862 }