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