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