src/gromacs/fileio/tngio.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014,2015,2018, 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 "gmxpre.h"
  36
  37 #include "tngio.h"
  38
  39 #include "config.h"
  40
  41 #include <cmath>
  42
  43 #include <algorithm>
  44 #include <iterator>
  45 #include <memory>
  46 #include <vector>
  47
  48 #if GMX_USE_TNG
  49 #include "tng/tng_io.h"
  50 #endif
  51
  52 #include "gromacs/math/units.h"
  53 #include "gromacs/math/utilities.h"
  54 #include "gromacs/mdtypes/inputrec.h"
  55 #include "gromacs/topology/ifunc.h"
  56 #include "gromacs/topology/topology.h"
  57 #include "gromacs/trajectory/trajectoryframe.h"
  58 #include "gromacs/utility/basedefinitions.h"
  59 #include "gromacs/utility/baseversion.h"
  60 #include "gromacs/utility/fatalerror.h"
  61 #include "gromacs/utility/futil.h"
  62 #include "gromacs/utility/gmxassert.h"
  63 #include "gromacs/utility/programcontext.h"
  64 #include "gromacs/utility/smalloc.h"
  65 #include "gromacs/utility/sysinfo.h"
  66 #include "gromacs/utility/unique_cptr.h"
  67
  68 #if !GMX_USE_TNG
  69 using tng_trajectory_t = void *;
  70 #endif
  71
  72 /*! \brief Gromacs Wrapper around tng datatype
  73  *
  74  * This could in principle hold any GROMACS-specific requirements not yet
  75  * implemented in or not relevant to the TNG library itself. However, for now
  76  * we only use it to handle some shortcomings we have discovered, where the TNG
  77  * API itself is a bit fragile and can end up overwriting data if called several
  78  * times with the same frame number.
  79  * The logic to determine the time per step was also a bit fragile. This is not
  80  * critical, but since we anyway need a wrapper for ensuring unique frame
  81  * numbers, we can also use it to store the time of the first step and use that
  82  * to derive a slightly better/safer estimate of the time per step.
  83  *
  84  * At some future point where we have a second-generation TNG API we should
  85  * consider removing this again.
  86  */
  87 struct gmx_tng_trajectory
  88 {
  89     tng_trajectory_t tng;                  //!< Actual TNG handle (pointer)
  90     bool             lastStepDataIsValid;  //!< True if lastStep has been set
  91     std::int64_t     lastStep;             //!< Index/step used for last frame
  92     bool             lastTimeDataIsValid;  //!< True if lastTime has been set
  93     double           lastTime;             //!< Time of last frame (TNG unit is seconds)
  94     bool             timePerFrameIsSet;    //!< True if we have set the time per frame
  95     int              boxOutputInterval;    //!< Number of steps between the output of box size
  96     int              lambdaOutputInterval; //!< Number of steps between the output of lambdas
  97 };
  98
  99 #if GMX_USE_TNG
 100 static const char *modeToVerb(char mode)
 101 {
 102     const char *p;
 103     switch (mode)
 104     {
 105         case 'r':
 106             p = "reading";
 107             break;
 108         case 'w':
 109             p = "writing";
 110             break;
 111         case 'a':
 112             p = "appending";
 113             break;
 114         default:
 115             gmx_fatal(FARGS, "Invalid file opening mode %c", mode);
 116     }
 117     return p;
 118 }
 119 #endif
 120
 121 void gmx_tng_open(const char           *filename,
 122                   char                  mode,
 123                   gmx_tng_trajectory_t *gmx_tng)
 124 {
 125 #if GMX_USE_TNG
 126     /* First check whether we have to make a backup,
 127      * only for writing, not for read or append.
 128      */
 129     if (mode == 'w')
 130     {
 131         make_backup(filename);
 132     }
 133
 134     *gmx_tng                        = new gmx_tng_trajectory;
 135     (*gmx_tng)->lastStepDataIsValid = false;
 136     (*gmx_tng)->lastTimeDataIsValid = false;
 137     (*gmx_tng)->timePerFrameIsSet   = false;
 138     tng_trajectory_t * tng = &(*gmx_tng)->tng;
 139
 140     /* tng must not be pointing at already allocated memory.
 141      * Memory will be allocated by tng_util_trajectory_open() and must
 142      * later on be freed by tng_util_trajectory_close(). */
 143     if (TNG_SUCCESS != tng_util_trajectory_open(filename, mode, tng))
 144     {
 145         /* TNG does return more than one degree of error, but there is
 146            no use case for GROMACS handling the non-fatal errors
 147            gracefully. */
 148         gmx_fatal(FARGS,
 149                   "File I/O error while opening %s for %s",
 150                   filename,
 151                   modeToVerb(mode));
 152     }
 153
 154     if (mode == 'w' || mode == 'a')
 155     {
 156         char hostname[256];
 157         gmx_gethostname(hostname, 256);
 158         if (mode == 'w')
 159         {
 160             tng_first_computer_name_set(*tng, hostname);
 161         }
 162         else
 163         {
 164             tng_last_computer_name_set(*tng, hostname);
 165         }
 166
 167         char        programInfo[256];
 168         const char *precisionString = "";
 169 #if GMX_DOUBLE
 170         precisionString = " (double precision)";
 171 #endif
 172         sprintf(programInfo, "%.100s %.128s%.24s",
 173                 gmx::getProgramContext().displayName(),
 174                 gmx_version(), precisionString);
 175         if (mode == 'w')
 176         {
 177             tng_first_program_name_set(*tng, programInfo);
 178         }
 179         else
 180         {
 181             tng_last_program_name_set(*tng, programInfo);
 182         }
 183
 184         char username[256];
 185         if (!gmx_getusername(username, 256))
 186         {
 187             if (mode == 'w')
 188             {
 189                 tng_first_user_name_set(*tng, username);
 190             }
 191             else
 192             {
 193                 tng_last_user_name_set(*tng, username);
 194                 tng_file_headers_write(*tng, TNG_USE_HASH);
 195             }
 196         }
 197     }
 198 #else
 199     gmx_file("GROMACS was compiled without TNG support, cannot handle this file type");
 200     GMX_UNUSED_VALUE(filename);
 201     GMX_UNUSED_VALUE(mode);
 202     GMX_UNUSED_VALUE(gmx_tng);
 203 #endif
 204 }
 205
 206 void gmx_tng_close(gmx_tng_trajectory_t *gmx_tng)
 207 {
 208     /* We have to check that tng is set because
 209      * tng_util_trajectory_close wants to return a NULL in it, and
 210      * gives a fatal error if it is NULL. */
 211 #if GMX_USE_TNG
 212     if (gmx_tng == nullptr || *gmx_tng == nullptr)
 213     {
 214         return;
 215     }
 216     tng_trajectory_t * tng = &(*gmx_tng)->tng;
 217
 218     if (tng)
 219     {
 220         tng_util_trajectory_close(tng);
 221     }
 222     delete *gmx_tng;
 223     *gmx_tng = nullptr;
 224
 225 #else
 226     GMX_UNUSED_VALUE(gmx_tng);
 227 #endif
 228 }
 229
 230 #if GMX_USE_TNG
 231 static void addTngMoleculeFromTopology(gmx_tng_trajectory_t gmx_tng,
 232                                        const char          *moleculeName,
 233                                        const t_atoms       *atoms,
 234                                        int64_t              numMolecules,
 235                                        tng_molecule_t      *tngMol)
 236 {
 237     tng_trajectory_t tng      = gmx_tng->tng;
 238     tng_chain_t      tngChain = nullptr;
 239     tng_residue_t    tngRes   = nullptr;
 240
 241     if (tng_molecule_add(tng, moleculeName, tngMol) != TNG_SUCCESS)
 242     {
 243         gmx_file("Cannot add molecule to TNG molecular system.");
 244     }
 245
 246     for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++)
 247     {
 248         const t_atom *at = &atoms->atom[atomIndex];
 249         /* FIXME: Currently the TNG API can only add atoms belonging to a
 250          * residue and chain. Wait for TNG 2.0*/
 251         if (atoms->nres > 0)
 252         {
 253             const t_resinfo *resInfo        = &atoms->resinfo[at->resind];
 254             char             chainName[2]   = {resInfo->chainid, 0};
 255             tng_atom_t       tngAtom        = nullptr;
 256             t_atom          *prevAtom;
 257
 258             if (atomIndex > 0)
 259             {
 260                 prevAtom = &atoms->atom[atomIndex - 1];
 261             }
 262             else
 263             {
 264                 prevAtom = nullptr;
 265             }
 266
 267             /* If this is the first atom or if the residue changed add the
 268              * residue to the TNG molecular system. */
 269             if (!prevAtom || resInfo != &atoms->resinfo[prevAtom->resind])
 270             {
 271                 /* If this is the first atom or if the chain changed add
 272                  * the chain to the TNG molecular system. */
 273                 if (!prevAtom || resInfo->chainid !=
 274                     atoms->resinfo[prevAtom->resind].chainid)
 275                 {
 276                     tng_molecule_chain_add(tng, *tngMol, chainName,
 277                                            &tngChain);
 278                 }
 279                 /* FIXME: When TNG supports both residue index and residue
 280                  * number the latter should be used. Wait for TNG 2.0*/
 281                 tng_chain_residue_add(tng, tngChain, *resInfo->name, &tngRes);
 282             }
 283             tng_residue_atom_add(tng, tngRes, *(atoms->atomname[atomIndex]), *(atoms->atomtype[atomIndex]), &tngAtom);
 284         }
 285     }
 286     tng_molecule_cnt_set(tng, *tngMol, numMolecules);
 287 }
 288
 289 void gmx_tng_add_mtop(gmx_tng_trajectory_t  gmx_tng,
 290                       const gmx_mtop_t     *mtop)
 291 {
 292     int                i;
 293     int                j;
 294     std::vector<real>  atomCharges;
 295     std::vector<real>  atomMasses;
 296     const t_ilist     *ilist;
 297     tng_bond_t         tngBond;
 298     char               datatype;
 299
 300     tng_trajectory_t   tng = gmx_tng->tng;
 301
 302     if (!mtop)
 303     {
 304         /* No topology information available to add. */
 305         return;
 306     }
 307
 308 #if GMX_DOUBLE
 309     datatype = TNG_DOUBLE_DATA;
 310 #else
 311     datatype = TNG_FLOAT_DATA;
 312 #endif
 313
 314     atomCharges.reserve(mtop->natoms);
 315     atomMasses.reserve(mtop->natoms);
 316
 317     for (const gmx_molblock_t &molBlock :  mtop->molblock)
 318     {
 319         tng_molecule_t       tngMol  = nullptr;
 320         const gmx_moltype_t *molType = &mtop->moltype[molBlock.type];
 321
 322         /* Add a molecule to the TNG trajectory with the same name as the
 323          * current molecule. */
 324         addTngMoleculeFromTopology(gmx_tng,
 325                                    *(molType->name),
 326                                    &molType->atoms,
 327                                    molBlock.nmol,
 328                                    &tngMol);
 329
 330         /* Bonds have to be deduced from interactions (constraints etc). Different
 331          * interactions have different sets of parameters. */
 332         /* Constraints are specified using two atoms */
 333         for (i = 0; i < F_NRE; i++)
 334         {
 335             if (IS_CHEMBOND(i))
 336             {
 337                 ilist = &molType->ilist[i];
 338                 if (ilist)
 339                 {
 340                     j = 1;
 341                     while (j < ilist->nr)
 342                     {
 343                         tng_molecule_bond_add(tng, tngMol, ilist->iatoms[j], ilist->iatoms[j+1], &tngBond);
 344                         j += 3;
 345                     }
 346                 }
 347             }
 348         }
 349         /* Settle is described using three atoms */
 350         ilist = &molType->ilist[F_SETTLE];
 351         if (ilist)
 352         {
 353             j = 1;
 354             while (j < ilist->nr)
 355             {
 356                 tng_molecule_bond_add(tng, tngMol, ilist->iatoms[j], ilist->iatoms[j+1], &tngBond);
 357                 tng_molecule_bond_add(tng, tngMol, ilist->iatoms[j], ilist->iatoms[j+2], &tngBond);
 358                 j += 4;
 359             }
 360         }
 361         /* First copy atom charges and masses, first atom by atom and then copy the memory for the molecule instances.
 362          * FIXME: Atom B state data should also be written to TNG (v 2.0?) */
 363         for (int atomCounter = 0; atomCounter < molType->atoms.nr; atomCounter++)
 364         {
 365             atomCharges.push_back(molType->atoms.atom[atomCounter].q);
 366             atomMasses.push_back(molType->atoms.atom[atomCounter].m);
 367         }
 368         for (int molCounter = 1; molCounter < molBlock.nmol; molCounter++)
 369         {
 370             std::copy_n(atomCharges.end() - molType->atoms.nr, molType->atoms.nr, std::back_inserter(atomCharges));
 371             std::copy_n(atomMasses.end()  - molType->atoms.nr, molType->atoms.nr, std::back_inserter(atomMasses));
 372         }
 373     }
 374     /* Write the TNG data blocks. */
 375     tng_particle_data_block_add(tng, TNG_TRAJ_PARTIAL_CHARGES, "PARTIAL CHARGES",
 376                                 datatype, TNG_NON_TRAJECTORY_BLOCK,
 377                                 1, 1, 1, 0, mtop->natoms,
 378                                 TNG_GZIP_COMPRESSION, atomCharges.data());
 379     tng_particle_data_block_add(tng, TNG_TRAJ_MASSES, "ATOM MASSES",
 380                                 datatype, TNG_NON_TRAJECTORY_BLOCK,
 381                                 1, 1, 1, 0, mtop->natoms,
 382                                 TNG_GZIP_COMPRESSION, atomMasses.data());
 383 }
 384
 385 /*! \libinternal \brief Compute greatest common divisor of n1 and n2
 386  * if they are positive.
 387  *
 388  * If only one of n1 and n2 is positive, then return it.
 389  * If neither n1 or n2 is positive, then return -1. */
 390 static int
 391 greatest_common_divisor_if_positive(int n1, int n2)
 392 {
 393     if (0 >= n1)
 394     {
 395         return (0 >= n2) ? -1 : n2;
 396     }
 397     if (0 >= n2)
 398     {
 399         return n1;
 400     }
 401
 402     /* We have a non-trivial greatest common divisor to compute. */
 403     return gmx_greatest_common_divisor(n1, n2);
 404 }
 405
 406 /* By default try to write 100 frames (of actual output) in each frame set.
 407  * This number is the number of outputs of the most frequently written data
 408  * type per frame set.
 409  * TODO for 5.1: Verify that 100 frames per frame set is efficient for most
 410  * setups regarding compression efficiency and compression time. Make this
 411  * a hidden command-line option? */
 412 const int defaultFramesPerFrameSet = 100;
 413
 414 /*! \libinternal \brief  Set the number of frames per frame
 415  * set according to output intervals.
 416  * The default is that 100 frames are written of the data
 417  * that is written most often. */
 418 static void tng_set_frames_per_frame_set(gmx_tng_trajectory_t  gmx_tng,
 419                                          const gmx_bool        bUseLossyCompression,
 420                                          const t_inputrec     *ir)
 421 {
 422     int              gcd = -1;
 423     tng_trajectory_t tng = gmx_tng->tng;
 424
 425     /* Set the number of frames per frame set to contain at least
 426      * defaultFramesPerFrameSet of the lowest common denominator of
 427      * the writing interval of positions and velocities. */
 428     /* FIXME after 5.0: consider nstenergy also? */
 429     if (bUseLossyCompression)
 430     {
 431         gcd = ir->nstxout_compressed;
 432     }
 433     else
 434     {
 435         gcd = greatest_common_divisor_if_positive(ir->nstxout, ir->nstvout);
 436         gcd = greatest_common_divisor_if_positive(gcd, ir->nstfout);
 437     }
 438     if (0 >= gcd)
 439     {
 440         return;
 441     }
 442
 443     tng_num_frames_per_frame_set_set(tng, gcd * defaultFramesPerFrameSet);
 444 }
 445
 446 /*! \libinternal \brief Set the data-writing intervals, and number of
 447  * frames per frame set */
 448 static void set_writing_intervals(gmx_tng_trajectory_t  gmx_tng,
 449                                   const gmx_bool        bUseLossyCompression,
 450                                   const t_inputrec     *ir)
 451 {
 452     tng_trajectory_t tng = gmx_tng->tng;
 453
 454     /* Define pointers to specific writing functions depending on if we
 455      * write float or double data */
 456     typedef tng_function_status (*set_writing_interval_func_pointer)(tng_trajectory_t,
 457                                                                      const int64_t,
 458                                                                      const int64_t,
 459                                                                      const int64_t,
 460                                                                      const char*,
 461                                                                      const char,
 462                                                                      const char);
 463 #if GMX_DOUBLE
 464     set_writing_interval_func_pointer set_writing_interval = tng_util_generic_write_interval_double_set;
 465 #else
 466     set_writing_interval_func_pointer set_writing_interval = tng_util_generic_write_interval_set;
 467 #endif
 468     int  xout, vout, fout;
 469     int  gcd = -1, lowest = -1;
 470     char compression;
 471
 472     tng_set_frames_per_frame_set(gmx_tng, bUseLossyCompression, ir);
 473
 474     if (bUseLossyCompression)
 475     {
 476         xout        = ir->nstxout_compressed;
 477
 478         /* If there is no uncompressed coordinate output write forces
 479            and velocities to the compressed tng file. */
 480         if (ir->nstxout)
 481         {
 482             vout        = 0;
 483             fout        = 0;
 484         }
 485         else
 486         {
 487             vout        = ir->nstvout;
 488             fout        = ir->nstfout;
 489         }
 490         compression = TNG_TNG_COMPRESSION;
 491     }
 492     else
 493     {
 494         xout        = ir->nstxout;
 495         vout        = ir->nstvout;
 496         fout        = ir->nstfout;
 497         compression = TNG_GZIP_COMPRESSION;
 498     }
 499     if (xout)
 500     {
 501         set_writing_interval(tng, xout, 3, TNG_TRAJ_POSITIONS,
 502                              "POSITIONS", TNG_PARTICLE_BLOCK_DATA,
 503                              compression);
 504         /* TODO: if/when we write energies to TNG also, reconsider how
 505          * and when box information is written, because GROMACS
 506          * behaviour pre-5.0 was to write the box with every
 507          * trajectory frame and every energy frame, and probably
 508          * people depend on this. */
 509
 510         gcd = greatest_common_divisor_if_positive(gcd, xout);
 511         if (lowest < 0 || xout < lowest)
 512         {
 513             lowest = xout;
 514         }
 515     }
 516     if (vout)
 517     {
 518         set_writing_interval(tng, vout, 3, TNG_TRAJ_VELOCITIES,
 519                              "VELOCITIES", TNG_PARTICLE_BLOCK_DATA,
 520                              compression);
 521
 522         gcd = greatest_common_divisor_if_positive(gcd, vout);
 523         if (lowest < 0 || vout < lowest)
 524         {
 525             lowest = vout;
 526         }
 527     }
 528     if (fout)
 529     {
 530         set_writing_interval(tng, fout, 3, TNG_TRAJ_FORCES,
 531                              "FORCES", TNG_PARTICLE_BLOCK_DATA,
 532                              TNG_GZIP_COMPRESSION);
 533
 534         gcd = greatest_common_divisor_if_positive(gcd, fout);
 535         if (lowest < 0 || fout < lowest)
 536         {
 537             lowest = fout;
 538         }
 539     }
 540     if (gcd > 0)
 541     {
 542         /* Lambdas and box shape written at an interval of the lowest common
 543            denominator of other output */
 544         set_writing_interval(tng, gcd, 1, TNG_GMX_LAMBDA,
 545                              "LAMBDAS", TNG_NON_PARTICLE_BLOCK_DATA,
 546                              TNG_GZIP_COMPRESSION);
 547
 548         set_writing_interval(tng, gcd, 9, TNG_TRAJ_BOX_SHAPE,
 549                              "BOX SHAPE", TNG_NON_PARTICLE_BLOCK_DATA,
 550                              TNG_GZIP_COMPRESSION);
 551         gmx_tng->lambdaOutputInterval = gcd;
 552         gmx_tng->boxOutputInterval    = gcd;
 553         if (gcd < lowest / 10)
 554         {
 555             gmx_warning("The lowest common denominator of trajectory output is "
 556                         "every %d step(s), whereas the shortest output interval "
 557                         "is every %d steps.", gcd, lowest);
 558         }
 559     }
 560 }
 561 #endif
 562
 563 void gmx_tng_prepare_md_writing(gmx_tng_trajectory_t  gmx_tng,
 564                                 const gmx_mtop_t     *mtop,
 565                                 const t_inputrec     *ir)
 566 {
 567 #if GMX_USE_TNG
 568     gmx_tng_add_mtop(gmx_tng, mtop);
 569     set_writing_intervals(gmx_tng, FALSE, ir);
 570     tng_time_per_frame_set(gmx_tng->tng, ir->delta_t * PICO);
 571     gmx_tng->timePerFrameIsSet = true;
 572 #else
 573     GMX_UNUSED_VALUE(gmx_tng);
 574     GMX_UNUSED_VALUE(mtop);
 575     GMX_UNUSED_VALUE(ir);
 576 #endif
 577 }
 578
 579 #if GMX_USE_TNG
 580 /* Check if all atoms in the molecule system are selected
 581  * by a selection group of type specified by gtype. */
 582 static gmx_bool all_atoms_selected(const gmx_mtop_t *mtop,
 583                                    const int         gtype)
 584 {
 585     /* Iterate through all atoms in the molecule system and
 586      * check if they belong to a selection group of the
 587      * requested type. */
 588     int i = 0;
 589     for (const gmx_molblock_t &molBlock : mtop->molblock)
 590     {
 591         const gmx_moltype_t &molType = mtop->moltype[molBlock.type];
 592         const t_atoms       &atoms   = molType.atoms;
 593
 594         for (int j = 0; j < molBlock.nmol; j++)
 595         {
 596             for (int atomIndex = 0; atomIndex < atoms.nr; atomIndex++, i++)
 597             {
 598                 if (ggrpnr(&mtop->groups, gtype, i) != 0)
 599                 {
 600                     return FALSE;
 601                 }
 602             }
 603         }
 604     }
 605     return TRUE;
 606 }
 607
 608 /* Create TNG molecules which will represent each of the selection
 609  * groups for writing. But do that only if there is actually a
 610  * specified selection group and it is not the whole system.
 611  * TODO: Currently the only selection that is taken into account
 612  * is egcCompressedX, but other selections should be added when
 613  * e.g. writing energies is implemented.
 614  */
 615 static void add_selection_groups(gmx_tng_trajectory_t  gmx_tng,
 616                                  const gmx_mtop_t     *mtop)
 617 {
 618     const t_atoms           *atoms;
 619     const t_atom            *at;
 620     const t_resinfo         *resInfo;
 621     const t_ilist           *ilist;
 622     int                      nameIndex;
 623     int                      atom_offset = 0;
 624     tng_molecule_t           mol, iterMol;
 625     tng_chain_t              chain;
 626     tng_residue_t            res;
 627     tng_atom_t               atom;
 628     tng_bond_t               tngBond;
 629     int64_t                  nMols;
 630     char                    *groupName;
 631     tng_trajectory_t         tng = gmx_tng->tng;
 632
 633     /* TODO: When the TNG molecules block is more flexible TNG selection
 634      * groups should not need all atoms specified. It should be possible
 635      * just to specify what molecules are selected (and/or which atoms in
 636      * the molecule) and how many (if applicable). */
 637
 638     /* If no atoms are selected we do not need to create a
 639      * TNG selection group molecule. */
 640     if (mtop->groups.ngrpnr[egcCompressedX] == 0)
 641     {
 642         return;
 643     }
 644
 645     /* If all atoms are selected we do not have to create a selection
 646      * group molecule in the TNG molecule system. */
 647     if (all_atoms_selected(mtop, egcCompressedX))
 648     {
 649         return;
 650     }
 651
 652     /* The name of the TNG molecule containing the selection group is the
 653      * same as the name of the selection group. */
 654     nameIndex = *mtop->groups.grps[egcCompressedX].nm_ind;
 655     groupName = *mtop->groups.grpname[nameIndex];
 656
 657     tng_molecule_alloc(tng, &mol);
 658     tng_molecule_name_set(tng, mol, groupName);
 659     tng_molecule_chain_add(tng, mol, "", &chain);
 660     int i = 0;
 661     for (const gmx_molblock_t &molBlock :  mtop->molblock)
 662     {
 663         const gmx_moltype_t &molType = mtop->moltype[molBlock.type];
 664
 665         atoms = &molType.atoms;
 666
 667         for (int j = 0; j < molBlock.nmol; j++)
 668         {
 669             bool bAtomsAdded = FALSE;
 670             for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++, i++)
 671             {
 672                 char *res_name;
 673                 int   res_id;
 674
 675                 if (ggrpnr(&mtop->groups, egcCompressedX, i) != 0)
 676                 {
 677                     continue;
 678                 }
 679                 at = &atoms->atom[atomIndex];
 680                 if (atoms->nres > 0)
 681                 {
 682                     resInfo = &atoms->resinfo[at->resind];
 683                     /* FIXME: When TNG supports both residue index and residue
 684                      * number the latter should be used. */
 685                     res_name = *resInfo->name;
 686                     res_id   = at->resind + 1;
 687                 }
 688                 else
 689                 {
 690                     res_name = const_cast<char*>("");
 691                     res_id   = 0;
 692                 }
 693                 if (tng_chain_residue_find(tng, chain, res_name, res_id, &res)
 694                     != TNG_SUCCESS)
 695                 {
 696                     /* Since there is ONE chain for selection groups do not keep the
 697                      * original residue IDs - otherwise there might be conflicts. */
 698                     tng_chain_residue_add(tng, chain, res_name, &res);
 699                 }
 700                 tng_residue_atom_w_id_add(tng, res, *(atoms->atomname[atomIndex]),
 701                                           *(atoms->atomtype[atomIndex]),
 702                                           atom_offset + atomIndex, &atom);
 703                 bAtomsAdded = TRUE;
 704             }
 705             /* Add bonds. */
 706             if (bAtomsAdded)
 707             {
 708                 for (int k = 0; k < F_NRE; k++)
 709                 {
 710                     if (IS_CHEMBOND(k))
 711                     {
 712                         ilist = &molType.ilist[k];
 713                         if (ilist)
 714                         {
 715                             int l = 1;
 716                             while (l < ilist->nr)
 717                             {
 718                                 int atom1, atom2;
 719                                 atom1 = ilist->iatoms[l] + atom_offset;
 720                                 atom2 = ilist->iatoms[l+1] + atom_offset;
 721                                 if (ggrpnr(&mtop->groups, egcCompressedX, atom1) == 0 &&
 722                                     ggrpnr(&mtop->groups, egcCompressedX, atom2) == 0)
 723                                 {
 724                                     tng_molecule_bond_add(tng, mol, ilist->iatoms[l],
 725                                                           ilist->iatoms[l+1], &tngBond);
 726                                 }
 727                                 l += 3;
 728                             }
 729                         }
 730                     }
 731                 }
 732                 /* Settle is described using three atoms */
 733                 ilist = &molType.ilist[F_SETTLE];
 734                 if (ilist)
 735                 {
 736                     int l = 1;
 737                     while (l < ilist->nr)
 738                     {
 739                         int atom1, atom2, atom3;
 740                         atom1 = ilist->iatoms[l] + atom_offset;
 741                         atom2 = ilist->iatoms[l+1] + atom_offset;
 742                         atom3 = ilist->iatoms[l+2] + atom_offset;
 743                         if (ggrpnr(&mtop->groups, egcCompressedX, atom1) == 0)
 744                         {
 745                             if (ggrpnr(&mtop->groups, egcCompressedX, atom2) == 0)
 746                             {
 747                                 tng_molecule_bond_add(tng, mol, atom1,
 748                                                       atom2, &tngBond);
 749                             }
 750                             if (ggrpnr(&mtop->groups, egcCompressedX, atom3) == 0)
 751                             {
 752                                 tng_molecule_bond_add(tng, mol, atom1,
 753                                                       atom3, &tngBond);
 754                             }
 755                         }
 756                         l += 4;
 757                     }
 758                 }
 759             }
 760             atom_offset += atoms->nr;
 761         }
 762     }
 763     tng_molecule_existing_add(tng, &mol);
 764     tng_molecule_cnt_set(tng, mol, 1);
 765     tng_num_molecule_types_get(tng, &nMols);
 766     for (int64_t k = 0; k < nMols; k++)
 767     {
 768         tng_molecule_of_index_get(tng, k, &iterMol);
 769         if (iterMol == mol)
 770         {
 771             continue;
 772         }
 773         tng_molecule_cnt_set(tng, iterMol, 0);
 774     }
 775 }
 776 #endif
 777
 778 void gmx_tng_set_compression_precision(gmx_tng_trajectory_t gmx_tng,
 779                                        real                 prec)
 780 {
 781 #if GMX_USE_TNG
 782     tng_compression_precision_set(gmx_tng->tng, prec);
 783 #else
 784     GMX_UNUSED_VALUE(gmx_tng);
 785     GMX_UNUSED_VALUE(prec);
 786 #endif
 787 }
 788
 789 void gmx_tng_prepare_low_prec_writing(gmx_tng_trajectory_t  gmx_tng,
 790                                       const gmx_mtop_t     *mtop,
 791                                       const t_inputrec     *ir)
 792 {
 793 #if GMX_USE_TNG
 794     gmx_tng_add_mtop(gmx_tng, mtop);
 795     add_selection_groups(gmx_tng, mtop);
 796     set_writing_intervals(gmx_tng, TRUE, ir);
 797     tng_time_per_frame_set(gmx_tng->tng, ir->delta_t * PICO);
 798     gmx_tng->timePerFrameIsSet = true;
 799     gmx_tng_set_compression_precision(gmx_tng, ir->x_compression_precision);
 800 #else
 801     GMX_UNUSED_VALUE(gmx_tng);
 802     GMX_UNUSED_VALUE(mtop);
 803     GMX_UNUSED_VALUE(ir);
 804 #endif
 805 }
 806
 807 void gmx_fwrite_tng(gmx_tng_trajectory_t gmx_tng,
 808                     const gmx_bool       bUseLossyCompression,
 809                     int64_t              step,
 810                     real                 elapsedPicoSeconds,
 811                     real                 lambda,
 812                     const rvec          *box,
 813                     int                  nAtoms,
 814                     const rvec          *x,
 815                     const rvec          *v,
 816                     const rvec          *f)
 817 {
 818 #if GMX_USE_TNG
 819     typedef tng_function_status (*write_data_func_pointer)(tng_trajectory_t,
 820                                                            const int64_t,
 821                                                            const double,
 822                                                            const real*,
 823                                                            const int64_t,
 824                                                            const int64_t,
 825                                                            const char*,
 826                                                            const char,
 827                                                            const char);
 828 #if GMX_DOUBLE
 829     static write_data_func_pointer           write_data           = tng_util_generic_with_time_double_write;
 830 #else
 831     static write_data_func_pointer           write_data           = tng_util_generic_with_time_write;
 832 #endif
 833     double                                   elapsedSeconds = elapsedPicoSeconds * PICO;
 834     int64_t                                  nParticles;
 835     char                                     compression;
 836
 837
 838     if (!gmx_tng)
 839     {
 840         /* This function might get called when the type of the
 841            compressed trajectory is actually XTC. So we exit and move
 842            on. */
 843         return;
 844     }
 845     tng_trajectory_t tng = gmx_tng->tng;
 846
 847     // While the GROMACS interface to this routine specifies 'step', TNG itself
 848     // only uses 'frame index' internally, although it suggests that it's a good
 849     // idea to let this represent the step rather than just frame numbers.
 850     //
 851     // However, the way the GROMACS interface works, there are cases where
 852     // the step is simply not set, so all frames rather have step=0.
 853     // If we call the lower-level TNG routines with the same frame index
 854     // (which is set from the step) multiple times, things will go very
 855     // wrong (overwritten data).
 856     //
 857     // To avoid this, we require the step value to always be larger than the
 858     // previous value, and if this is not true we simply set it to a value
 859     // one unit larger than the previous step.
 860     //
 861     // Note: We do allow the user to specify any crazy value the want for the
 862     // first step. If it is "not set", GROMACS will have used the value 0.
 863     if (gmx_tng->lastStepDataIsValid && step <= gmx_tng->lastStep)
 864     {
 865         step = gmx_tng->lastStep + 1;
 866     }
 867
 868     // Now that we have fixed the potentially incorrect step, we can also set
 869     // the time per frame if it was not already done, and we have
 870     // step/time information from the last frame so it is possible to calculate it.
 871     // Note that we do allow the time to be the same for all steps, or even
 872     // decreasing. In that case we will simply say the time per step is 0
 873     // or negative. A bit strange, but a correct representation of the data :-)
 874     if (!gmx_tng->timePerFrameIsSet && gmx_tng->lastTimeDataIsValid && gmx_tng->lastStepDataIsValid)
 875     {
 876         double       deltaTime = elapsedSeconds - gmx_tng->lastTime;
 877         std::int64_t deltaStep = step - gmx_tng->lastStep;
 878         tng_time_per_frame_set(tng, deltaTime / deltaStep );
 879         gmx_tng->timePerFrameIsSet = true;
 880     }
 881
 882     tng_num_particles_get(tng, &nParticles);
 883     if (nAtoms != static_cast<int>(nParticles))
 884     {
 885         tng_implicit_num_particles_set(tng, nAtoms);
 886     }
 887
 888     if (bUseLossyCompression)
 889     {
 890         compression = TNG_TNG_COMPRESSION;
 891     }
 892     else
 893     {
 894         compression = TNG_GZIP_COMPRESSION;
 895     }
 896
 897     /* The writing is done using write_data, which writes float or double
 898      * depending on the GROMACS compilation. */
 899     if (x)
 900     {
 901         GMX_ASSERT(box, "Need a non-NULL box if positions are written");
 902
 903         if (write_data(tng, step, elapsedSeconds,
 904                        reinterpret_cast<const real *>(x),
 905                        3, TNG_TRAJ_POSITIONS, "POSITIONS",
 906                        TNG_PARTICLE_BLOCK_DATA,
 907                        compression) != TNG_SUCCESS)
 908         {
 909             gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 910         }
 911     }
 912
 913     if (v)
 914     {
 915         if (write_data(tng, step, elapsedSeconds,
 916                        reinterpret_cast<const real *>(v),
 917                        3, TNG_TRAJ_VELOCITIES, "VELOCITIES",
 918                        TNG_PARTICLE_BLOCK_DATA,
 919                        compression) != TNG_SUCCESS)
 920         {
 921             gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 922         }
 923     }
 924
 925     if (f)
 926     {
 927         /* TNG-MF1 compression only compresses positions and velocities. Use lossless
 928          * compression for forces regardless of output mode */
 929         if (write_data(tng, step, elapsedSeconds,
 930                        reinterpret_cast<const real *>(f),
 931                        3, TNG_TRAJ_FORCES, "FORCES",
 932                        TNG_PARTICLE_BLOCK_DATA,
 933                        TNG_GZIP_COMPRESSION) != TNG_SUCCESS)
 934         {
 935             gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 936         }
 937     }
 938
 939     if (box)
 940     {
 941         /* TNG-MF1 compression only compresses positions and velocities. Use lossless
 942          * compression for box shape regardless of output mode */
 943         if (write_data(tng, step, elapsedSeconds,
 944                        reinterpret_cast<const real *>(box),
 945                        9, TNG_TRAJ_BOX_SHAPE, "BOX SHAPE",
 946                        TNG_NON_PARTICLE_BLOCK_DATA,
 947                        TNG_GZIP_COMPRESSION) != TNG_SUCCESS)
 948         {
 949             gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 950         }
 951     }
 952
 953     if (lambda >= 0)
 954     {
 955         /* TNG-MF1 compression only compresses positions and velocities. Use lossless
 956          * compression for lambda regardless of output mode */
 957         if (write_data(tng, step, elapsedSeconds,
 958                        reinterpret_cast<const real *>(&lambda),
 959                        1, TNG_GMX_LAMBDA, "LAMBDAS",
 960                        TNG_NON_PARTICLE_BLOCK_DATA,
 961                        TNG_GZIP_COMPRESSION) != TNG_SUCCESS)
 962         {
 963             gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 964         }
 965     }
 966
 967     // Update the data in the wrapper
 968
 969     gmx_tng->lastStepDataIsValid = true;
 970     gmx_tng->lastStep            = step;
 971     gmx_tng->lastTimeDataIsValid = true;
 972     gmx_tng->lastTime            = elapsedSeconds;
 973 #else
 974     GMX_UNUSED_VALUE(gmx_tng);
 975     GMX_UNUSED_VALUE(bUseLossyCompression);
 976     GMX_UNUSED_VALUE(step);
 977     GMX_UNUSED_VALUE(elapsedPicoSeconds);
 978     GMX_UNUSED_VALUE(lambda);
 979     GMX_UNUSED_VALUE(box);
 980     GMX_UNUSED_VALUE(nAtoms);
 981     GMX_UNUSED_VALUE(x);
 982     GMX_UNUSED_VALUE(v);
 983     GMX_UNUSED_VALUE(f);
 984 #endif
 985 }
 986
 987 void fflush_tng(gmx_tng_trajectory_t gmx_tng)
 988 {
 989 #if GMX_USE_TNG
 990     if (!gmx_tng)
 991     {
 992         return;
 993     }
 994     tng_frame_set_premature_write(gmx_tng->tng, TNG_USE_HASH);
 995 #else
 996     GMX_UNUSED_VALUE(gmx_tng);
 997 #endif
 998 }
 999
1000 float gmx_tng_get_time_of_final_frame(gmx_tng_trajectory_t gmx_tng)
1001 {
1002 #if GMX_USE_TNG
1003     int64_t          nFrames;
1004     double           time;
1005     float            fTime;
1006     tng_trajectory_t tng = gmx_tng->tng;
1007
1008     tng_num_frames_get(tng, &nFrames);
1009     tng_util_time_of_frame_get(tng, nFrames - 1, &time);
1010
1011     fTime = time / PICO;
1012     return fTime;
1013 #else
1014     GMX_UNUSED_VALUE(gmx_tng);
1015     return -1.0;
1016 #endif
1017 }
1018
1019 void gmx_prepare_tng_writing(const char              *filename,
1020                              char                     mode,
1021                              gmx_tng_trajectory_t    *gmx_tng_input,
1022                              gmx_tng_trajectory_t    *gmx_tng_output,
1023                              int                      nAtoms,
1024                              const gmx_mtop_t        *mtop,
1025                              const int               *index,
1026                              const char              *indexGroupName)
1027 {
1028 #if GMX_USE_TNG
1029     tng_trajectory_t   *input  = (gmx_tng_input && *gmx_tng_input) ? &(*gmx_tng_input)->tng : nullptr;
1030     /* FIXME after 5.0: Currently only standard block types are read */
1031     const int           defaultNumIds              = 5;
1032     static int64_t      fallbackIds[defaultNumIds] =
1033     {
1034         TNG_TRAJ_BOX_SHAPE, TNG_TRAJ_POSITIONS,
1035         TNG_TRAJ_VELOCITIES, TNG_TRAJ_FORCES,
1036         TNG_GMX_LAMBDA
1037     };
1038     static char         fallbackNames[defaultNumIds][32] =
1039     {
1040         "BOX SHAPE", "POSITIONS", "VELOCITIES",
1041         "FORCES", "LAMBDAS"
1042     };
1043
1044     typedef tng_function_status (*set_writing_interval_func_pointer)(tng_trajectory_t,
1045                                                                      const int64_t,
1046                                                                      const int64_t,
1047                                                                      const int64_t,
1048                                                                      const char*,
1049                                                                      const char,
1050                                                                      const char);
1051 #if GMX_DOUBLE
1052     set_writing_interval_func_pointer set_writing_interval = tng_util_generic_write_interval_double_set;
1053 #else
1054     set_writing_interval_func_pointer set_writing_interval = tng_util_generic_write_interval_set;
1055 #endif
1056
1057     gmx_tng_open(filename, mode, gmx_tng_output);
1058     tng_trajectory_t *output = &(*gmx_tng_output)->tng;
1059
1060     /* Do we have an input file in TNG format? If so, then there's
1061        more data we can copy over, rather than having to improvise. */
1062     if (gmx_tng_input && *gmx_tng_input)
1063     {
1064         /* Set parameters (compression, time per frame, molecule
1065          * information, number of frames per frame set and writing
1066          * intervals of positions, box shape and lambdas) of the
1067          * output tng container based on their respective values int
1068          * the input tng container */
1069         double      time, compression_precision;
1070         int64_t     n_frames_per_frame_set, interval = -1;
1071
1072         tng_compression_precision_get(*input, &compression_precision);
1073         tng_compression_precision_set(*output, compression_precision);
1074         // TODO make this configurable in a future version
1075         char compression_type = TNG_TNG_COMPRESSION;
1076
1077         tng_molecule_system_copy(*input, *output);
1078
1079         tng_time_per_frame_get(*input, &time);
1080         tng_time_per_frame_set(*output, time);
1081         // Since we have copied the value from the input TNG we should not change it again
1082         (*gmx_tng_output)->timePerFrameIsSet = true;
1083
1084         tng_num_frames_per_frame_set_get(*input, &n_frames_per_frame_set);
1085         tng_num_frames_per_frame_set_set(*output, n_frames_per_frame_set);
1086
1087         for (int i = 0; i < defaultNumIds; i++)
1088         {
1089             if (tng_data_get_stride_length(*input, fallbackIds[i], -1, &interval)
1090                 == TNG_SUCCESS)
1091             {
1092                 switch (fallbackIds[i])
1093                 {
1094                     case TNG_TRAJ_POSITIONS:
1095                     case TNG_TRAJ_VELOCITIES:
1096                         set_writing_interval(*output, interval, 3, fallbackIds[i],
1097                                              fallbackNames[i], TNG_PARTICLE_BLOCK_DATA,
1098                                              compression_type);
1099                         break;
1100                     case TNG_TRAJ_FORCES:
1101                         set_writing_interval(*output, interval, 3, fallbackIds[i],
1102                                              fallbackNames[i], TNG_PARTICLE_BLOCK_DATA,
1103                                              TNG_GZIP_COMPRESSION);
1104                         break;
1105                     case TNG_TRAJ_BOX_SHAPE:
1106                         set_writing_interval(*output, interval, 9, fallbackIds[i],
1107                                              fallbackNames[i], TNG_NON_PARTICLE_BLOCK_DATA,
1108                                              TNG_GZIP_COMPRESSION);
1109                         (*gmx_tng_output)->boxOutputInterval = interval;
1110                         break;
1111                     case TNG_GMX_LAMBDA:
1112                         set_writing_interval(*output, interval, 1, fallbackIds[i],
1113                                              fallbackNames[i], TNG_NON_PARTICLE_BLOCK_DATA,
1114                                              TNG_GZIP_COMPRESSION);
1115                         (*gmx_tng_output)->lambdaOutputInterval = interval;
1116                         break;
1117                     default:
1118                         continue;
1119                 }
1120             }
1121         }
1122
1123     }
1124     else
1125     {
1126         /* TODO after trjconv is modularized: fix this so the user can
1127            change precision when they are doing an operation where
1128            this makes sense, and not otherwise.
1129
1130            char compression = bUseLossyCompression ? TNG_TNG_COMPRESSION : TNG_GZIP_COMPRESSION;
1131            gmx_tng_set_compression_precision(*output, ndec2prec(nDecimalsOfPrecision));
1132          */
1133         gmx_tng_add_mtop(*gmx_tng_output, mtop);
1134         tng_num_frames_per_frame_set_set(*output, 1);
1135     }
1136
1137     if (index && nAtoms > 0)
1138     {
1139         gmx_tng_setup_atom_subgroup(*gmx_tng_output, nAtoms, index, indexGroupName);
1140     }
1141
1142     /* If for some reason there are more requested atoms than there are atoms in the
1143      * molecular system create a number of implicit atoms (without atom data) to
1144      * compensate for that. */
1145     if (nAtoms >= 0)
1146     {
1147         tng_implicit_num_particles_set(*output, nAtoms);
1148     }
1149 #else
1150     GMX_UNUSED_VALUE(filename);
1151     GMX_UNUSED_VALUE(mode);
1152     GMX_UNUSED_VALUE(gmx_tng_input);
1153     GMX_UNUSED_VALUE(gmx_tng_output);
1154     GMX_UNUSED_VALUE(nAtoms);
1155     GMX_UNUSED_VALUE(mtop);
1156     GMX_UNUSED_VALUE(index);
1157     GMX_UNUSED_VALUE(indexGroupName);
1158 #endif
1159 }
1160
1161 void gmx_write_tng_from_trxframe(gmx_tng_trajectory_t    gmx_tng_output,
1162                                  const t_trxframe       *frame,
1163                                  int                     natoms)
1164 {
1165 #if GMX_USE_TNG
1166     if (natoms < 0)
1167     {
1168         natoms = frame->natoms;
1169     }
1170     gmx_fwrite_tng(gmx_tng_output,
1171                    TRUE,
1172                    frame->step,
1173                    frame->time,
1174                    0,
1175                    frame->box,
1176                    natoms,
1177                    frame->x,
1178                    frame->v,
1179                    frame->f);
1180 #else
1181     GMX_UNUSED_VALUE(gmx_tng_output);
1182     GMX_UNUSED_VALUE(frame);
1183     GMX_UNUSED_VALUE(natoms);
1184 #endif
1185 }
1186
1187 namespace
1188 {
1189
1190 #if GMX_USE_TNG
1191 void
1192 convert_array_to_real_array(void       *from,
1193                             real       *to,
1194                             const float fact,
1195                             const int   nAtoms,
1196                             const int   nValues,
1197                             const char  datatype)
1198 {
1199     int        i, j;
1200
1201     const bool useDouble = GMX_DOUBLE;
1202     switch (datatype)
1203     {
1204         case TNG_FLOAT_DATA:
1205             if (!useDouble)
1206             {
1207                 if (fact == 1)
1208                 {
1209                     memcpy(to, from, nValues * sizeof(real) * nAtoms);
1210                 }
1211                 else
1212                 {
1213                     for (i = 0; i < nAtoms; i++)
1214                     {
1215                         for (j = 0; j < nValues; j++)
1216                         {
1217                             to[i*nValues+j] = reinterpret_cast<float *>(from)[i*nValues+j] * fact;
1218                         }
1219                     }
1220                 }
1221             }
1222             else
1223             {
1224                 for (i = 0; i < nAtoms; i++)
1225                 {
1226                     for (j = 0; j < nValues; j++)
1227                     {
1228                         to[i*nValues+j] = reinterpret_cast<float *>(from)[i*nValues+j] * fact;
1229                     }
1230                 }
1231             }
1232             break;
1233         case TNG_INT_DATA:
1234             for (i = 0; i < nAtoms; i++)
1235             {
1236                 for (j = 0; j < nValues; j++)
1237                 {
1238                     to[i*nValues+j] = reinterpret_cast<int64_t *>(from)[i*nValues+j] * fact;
1239                 }
1240             }
1241             break;
1242         case TNG_DOUBLE_DATA:
1243             if (sizeof(real) == sizeof(double))
1244             {
1245                 if (fact == 1)
1246                 {
1247                     memcpy(to, from, nValues * sizeof(real) * nAtoms);
1248                 }
1249                 else
1250                 {
1251                     for (i = 0; i < nAtoms; i++)
1252                     {
1253                         for (j = 0; j < nValues; j++)
1254                         {
1255                             to[i*nValues+j] = reinterpret_cast<double *>(from)[i*nValues+j] * fact;
1256                         }
1257                     }
1258                 }
1259             }
1260             else
1261             {
1262                 for (i = 0; i < nAtoms; i++)
1263                 {
1264                     for (j = 0; j < nValues; j++)
1265                     {
1266                         to[i*nValues+j] = reinterpret_cast<double *>(from)[i*nValues+j] * fact;
1267                     }
1268                 }
1269             }
1270             break;
1271         default:
1272             gmx_incons("Illegal datatype when converting values to a real array!");
1273     }
1274 }
1275
1276 real getDistanceScaleFactor(gmx_tng_trajectory_t in)
1277 {
1278     int64_t     exp = -1;
1279     real        distanceScaleFactor;
1280
1281     // TODO Hopefully, TNG 2.0 will do this kind of thing for us
1282     tng_distance_unit_exponential_get(in->tng, &exp);
1283
1284     // GROMACS expects distances in nm
1285     switch (exp)
1286     {
1287         case 9:
1288             distanceScaleFactor = NANO/NANO;
1289             break;
1290         case 10:
1291             distanceScaleFactor = NANO/ANGSTROM;
1292             break;
1293         default:
1294             distanceScaleFactor = pow(10.0, exp + 9.0);
1295     }
1296
1297     return distanceScaleFactor;
1298 }
1299 #endif
1300
1301 } // namespace
1302
1303 void gmx_tng_setup_atom_subgroup(gmx_tng_trajectory_t gmx_tng,
1304                                  const int            nind,
1305                                  const int           *ind,
1306                                  const char          *name)
1307 {
1308 #if GMX_USE_TNG
1309     int64_t                  nAtoms, cnt, nMols;
1310     tng_molecule_t           mol, iterMol;
1311     tng_chain_t              chain;
1312     tng_residue_t            res;
1313     tng_atom_t               atom;
1314     tng_function_status      stat;
1315     tng_trajectory_t         tng = gmx_tng->tng;
1316
1317     tng_num_particles_get(tng, &nAtoms);
1318
1319     if (nAtoms == nind)
1320     {
1321         return;
1322     }
1323
1324     stat = tng_molecule_find(tng, name, -1, &mol);
1325     if (stat == TNG_SUCCESS)
1326     {
1327         tng_molecule_num_atoms_get(tng, mol, &nAtoms);
1328         tng_molecule_cnt_get(tng, mol, &cnt);
1329         if (nAtoms == nind)
1330         {
1331             stat = TNG_SUCCESS;
1332         }
1333         else
1334         {
1335             stat = TNG_FAILURE;
1336         }
1337     }
1338     if (stat == TNG_FAILURE)
1339     {
1340         /* The indexed atoms are added to one separate molecule. */
1341         tng_molecule_alloc(tng, &mol);
1342         tng_molecule_name_set(tng, mol, name);
1343         tng_molecule_chain_add(tng, mol, "", &chain);
1344
1345         for (int i = 0; i < nind; i++)
1346         {
1347             char        temp_name[256], temp_type[256];
1348
1349             /* Try to retrieve the residue name of the atom */
1350             stat = tng_residue_name_of_particle_nr_get(tng, ind[i], temp_name, 256);
1351             if (stat != TNG_SUCCESS)
1352             {
1353                 temp_name[0] = '\0';
1354             }
1355             /* Check if the molecule of the selection already contains this residue */
1356             if (tng_chain_residue_find(tng, chain, temp_name, -1, &res)
1357                 != TNG_SUCCESS)
1358             {
1359                 tng_chain_residue_add(tng, chain, temp_name, &res);
1360             }
1361             /* Try to find the original name and type of the atom */
1362             stat = tng_atom_name_of_particle_nr_get(tng, ind[i], temp_name, 256);
1363             if (stat != TNG_SUCCESS)
1364             {
1365                 temp_name[0] = '\0';
1366             }
1367             stat = tng_atom_type_of_particle_nr_get(tng, ind[i], temp_type, 256);
1368             if (stat != TNG_SUCCESS)
1369             {
1370                 temp_type[0] = '\0';
1371             }
1372             tng_residue_atom_w_id_add(tng, res, temp_name, temp_type, ind[i], &atom);
1373         }
1374         tng_molecule_existing_add(tng, &mol);
1375     }
1376     /* Set the count of the molecule containing the selected atoms to 1 and all
1377      * other molecules to 0 */
1378     tng_molecule_cnt_set(tng, mol, 1);
1379     tng_num_molecule_types_get(tng, &nMols);
1380     for (int64_t k = 0; k < nMols; k++)
1381     {
1382         tng_molecule_of_index_get(tng, k, &iterMol);
1383         if (iterMol == mol)
1384         {
1385             continue;
1386         }
1387         tng_molecule_cnt_set(tng, iterMol, 0);
1388     }
1389 #else
1390     GMX_UNUSED_VALUE(gmx_tng);
1391     GMX_UNUSED_VALUE(nind);
1392     GMX_UNUSED_VALUE(ind);
1393     GMX_UNUSED_VALUE(name);
1394 #endif
1395 }
1396
1397 /* TODO: If/when TNG acquires the ability to copy data blocks without
1398  * uncompressing them, then this implemenation should be reconsidered.
1399  * Ideally, gmx trjconv -f a.tng -o b.tng -b 10 -e 20 would be fast
1400  * and lose no information. */
1401 gmx_bool gmx_read_next_tng_frame(gmx_tng_trajectory_t        gmx_tng_input,
1402                                  t_trxframe                 *fr,
1403                                  int64_t                    *requestedIds,
1404                                  int                         numRequestedIds)
1405 {
1406 #if GMX_USE_TNG
1407     tng_trajectory_t        input = gmx_tng_input->tng;
1408     gmx_bool                bOK   = TRUE;
1409     tng_function_status     stat;
1410     int64_t                 numberOfAtoms = -1, frameNumber = -1;
1411     int64_t                 nBlocks, blockId, *blockIds = nullptr, codecId;
1412     char                    datatype      = -1;
1413     void                   *values        = nullptr;
1414     double                  frameTime     = -1.0;
1415     int                     size, blockDependency;
1416     double                  prec;
1417     const int               defaultNumIds = 5;
1418     static int64_t          fallbackRequestedIds[defaultNumIds] =
1419     {
1420         TNG_TRAJ_BOX_SHAPE, TNG_TRAJ_POSITIONS,
1421         TNG_TRAJ_VELOCITIES, TNG_TRAJ_FORCES,
1422         TNG_GMX_LAMBDA
1423     };
1424
1425
1426     fr->bStep     = FALSE;
1427     fr->bTime     = FALSE;
1428     fr->bLambda   = FALSE;
1429     fr->bAtoms    = FALSE;
1430     fr->bPrec     = FALSE;
1431     fr->bX        = FALSE;
1432     fr->bV        = FALSE;
1433     fr->bF        = FALSE;
1434     fr->bBox      = FALSE;
1435
1436     /* If no specific IDs were requested read all block types that can
1437      * currently be interpreted */
1438     if (!requestedIds || numRequestedIds == 0)
1439     {
1440         numRequestedIds = defaultNumIds;
1441         requestedIds    = fallbackRequestedIds;
1442     }
1443
1444     stat = tng_num_particles_get(input, &numberOfAtoms);
1445     if (stat != TNG_SUCCESS)
1446     {
1447         gmx_file("Cannot determine number of atoms from TNG file.");
1448     }
1449     fr->natoms = numberOfAtoms;
1450
1451     bool nextFrameExists = gmx_get_tng_data_block_types_of_next_frame(gmx_tng_input,
1452                                                                       fr->step,
1453                                                                       numRequestedIds,
1454                                                                       requestedIds,
1455                                                                       &frameNumber,
1456                                                                       &nBlocks,
1457                                                                       &blockIds);
1458     gmx::unique_cptr<int64_t, gmx::free_wrapper> blockIdsGuard(blockIds);
1459     if (!nextFrameExists)
1460     {
1461         return FALSE;
1462     }
1463
1464     if (nBlocks == 0)
1465     {
1466         return FALSE;
1467     }
1468
1469     for (int64_t i = 0; i < nBlocks; i++)
1470     {
1471         blockId = blockIds[i];
1472         tng_data_block_dependency_get(input, blockId, &blockDependency);
1473         if (blockDependency & TNG_PARTICLE_DEPENDENT)
1474         {
1475             stat = tng_util_particle_data_next_frame_read(input,
1476                                                           blockId,
1477                                                           &values,
1478                                                           &datatype,
1479                                                           &frameNumber,
1480                                                           &frameTime);
1481         }
1482         else
1483         {
1484             stat = tng_util_non_particle_data_next_frame_read(input,
1485                                                               blockId,
1486                                                               &values,
1487                                                               &datatype,
1488                                                               &frameNumber,
1489                                                               &frameTime);
1490         }
1491         if (stat == TNG_CRITICAL)
1492         {
1493             gmx_file("Cannot read positions from TNG file.");
1494             return FALSE;
1495         }
1496         else if (stat == TNG_FAILURE)
1497         {
1498             continue;
1499         }
1500         switch (blockId)
1501         {
1502             case TNG_TRAJ_BOX_SHAPE:
1503                 switch (datatype)
1504                 {
1505                     case TNG_INT_DATA:
1506                         size = sizeof(int64_t);
1507                         break;
1508                     case TNG_FLOAT_DATA:
1509                         size = sizeof(float);
1510                         break;
1511                     case TNG_DOUBLE_DATA:
1512                         size = sizeof(double);
1513                         break;
1514                     default:
1515                         gmx_incons("Illegal datatype of box shape values!");
1516                 }
1517                 for (int i = 0; i < DIM; i++)
1518                 {
1519                     convert_array_to_real_array(reinterpret_cast<char *>(values) + size * i * DIM,
1520                                                 reinterpret_cast<real *>(fr->box[i]),
1521                                                 getDistanceScaleFactor(gmx_tng_input),
1522                                                 1,
1523                                                 DIM,
1524                                                 datatype);
1525                 }
1526                 fr->bBox = TRUE;
1527                 break;
1528             case TNG_TRAJ_POSITIONS:
1529                 srenew(fr->x, fr->natoms);
1530                 convert_array_to_real_array(values,
1531                                             reinterpret_cast<real *>(fr->x),
1532                                             getDistanceScaleFactor(gmx_tng_input),
1533                                             fr->natoms,
1534                                             DIM,
1535                                             datatype);
1536                 fr->bX = TRUE;
1537                 tng_util_frame_current_compression_get(input, blockId, &codecId, &prec);
1538                 /* This must be updated if/when more lossy compression methods are added */
1539                 if (codecId == TNG_TNG_COMPRESSION)
1540                 {
1541                     fr->prec  = prec;
1542                     fr->bPrec = TRUE;
1543                 }
1544                 break;
1545             case TNG_TRAJ_VELOCITIES:
1546                 srenew(fr->v, fr->natoms);
1547                 convert_array_to_real_array(values,
1548                                             reinterpret_cast<real *>(fr->v),
1549                                             getDistanceScaleFactor(gmx_tng_input),
1550                                             fr->natoms,
1551                                             DIM,
1552                                             datatype);
1553                 fr->bV = TRUE;
1554                 tng_util_frame_current_compression_get(input, blockId, &codecId, &prec);
1555                 /* This must be updated if/when more lossy compression methods are added */
1556                 if (codecId == TNG_TNG_COMPRESSION)
1557                 {
1558                     fr->prec  = prec;
1559                     fr->bPrec = TRUE;
1560                 }
1561                 break;
1562             case TNG_TRAJ_FORCES:
1563                 srenew(fr->f, fr->natoms);
1564                 convert_array_to_real_array(values,
1565                                             reinterpret_cast<real *>(fr->f),
1566                                             getDistanceScaleFactor(gmx_tng_input),
1567                                             fr->natoms,
1568                                             DIM,
1569                                             datatype);
1570                 fr->bF = TRUE;
1571                 break;
1572             case TNG_GMX_LAMBDA:
1573                 switch (datatype)
1574                 {
1575                     case TNG_FLOAT_DATA:
1576                         fr->lambda = *(reinterpret_cast<float *>(values));
1577                         break;
1578                     case TNG_DOUBLE_DATA:
1579                         fr->lambda = *(reinterpret_cast<double *>(values));
1580                         break;
1581                     default:
1582                         gmx_incons("Illegal datatype lambda value!");
1583                 }
1584                 fr->bLambda = TRUE;
1585                 break;
1586             default:
1587                 gmx_warning("Illegal block type! Currently GROMACS tools can only handle certain data types. Skipping block.");
1588         }
1589         /* values does not have to be freed before reading next frame. It will
1590          * be reallocated if it is not NULL. */
1591     }
1592
1593     fr->step  = frameNumber;
1594     fr->bStep = TRUE;
1595
1596     // Convert the time to ps
1597     fr->time  = frameTime / PICO;
1598     fr->bTime = (frameTime > 0);
1599
1600     // TODO This does not leak, but is not exception safe.
1601     /* values must be freed before leaving this function */
1602     sfree(values);
1603
1604     return bOK;
1605 #else
1606     GMX_UNUSED_VALUE(gmx_tng_input);
1607     GMX_UNUSED_VALUE(fr);
1608     GMX_UNUSED_VALUE(requestedIds);
1609     GMX_UNUSED_VALUE(numRequestedIds);
1610     return FALSE;
1611 #endif
1612 }
1613
1614 void gmx_print_tng_molecule_system(gmx_tng_trajectory_t gmx_tng_input,
1615                                    FILE                *stream)
1616 {
1617 #if GMX_USE_TNG
1618     int64_t             nMolecules, nChains, nResidues, nAtoms, nFramesRead;
1619     int64_t             strideLength, nParticlesRead, nValuesPerFrameRead, *molCntList;
1620     tng_molecule_t      molecule;
1621     tng_chain_t         chain;
1622     tng_residue_t       residue;
1623     tng_atom_t          atom;
1624     tng_function_status stat;
1625     char                str[256];
1626     char                varNAtoms;
1627     char                datatype;
1628     void               *data = nullptr;
1629     std::vector<real>   atomCharges;
1630     std::vector<real>   atomMasses;
1631     tng_trajectory_t    input = gmx_tng_input->tng;
1632
1633     tng_num_molecule_types_get(input, &nMolecules);
1634     tng_molecule_cnt_list_get(input, &molCntList);
1635     /* Can the number of particles change in the trajectory or is it constant? */
1636     tng_num_particles_variable_get(input, &varNAtoms);
1637
1638     for (int64_t i = 0; i < nMolecules; i++)
1639     {
1640         tng_molecule_of_index_get(input, i, &molecule);
1641         tng_molecule_name_get(input, molecule, str, 256);
1642         if (varNAtoms == TNG_CONSTANT_N_ATOMS)
1643         {
1644             if (static_cast<int>(molCntList[i]) == 0)
1645             {
1646                 continue;
1647             }
1648             fprintf(stream, "Molecule: %s, count: %d\n", str, static_cast<int>(molCntList[i]));
1649         }
1650         else
1651         {
1652             fprintf(stream, "Molecule: %s\n", str);
1653         }
1654         tng_molecule_num_chains_get(input, molecule, &nChains);
1655         if (nChains > 0)
1656         {
1657             for (int64_t j = 0; j < nChains; j++)
1658             {
1659                 tng_molecule_chain_of_index_get(input, molecule, j, &chain);
1660                 tng_chain_name_get(input, chain, str, 256);
1661                 fprintf(stream, "\tChain: %s\n", str);
1662                 tng_chain_num_residues_get(input, chain, &nResidues);
1663                 for (int64_t k = 0; k < nResidues; k++)
1664                 {
1665                     tng_chain_residue_of_index_get(input, chain, k, &residue);
1666                     tng_residue_name_get(input, residue, str, 256);
1667                     fprintf(stream, "\t\tResidue: %s\n", str);
1668                     tng_residue_num_atoms_get(input, residue, &nAtoms);
1669                     for (int64_t l = 0; l < nAtoms; l++)
1670                     {
1671                         tng_residue_atom_of_index_get(input, residue, l, &atom);
1672                         tng_atom_name_get(input, atom, str, 256);
1673                         fprintf(stream, "\t\t\tAtom: %s", str);
1674                         tng_atom_type_get(input, atom, str, 256);
1675                         fprintf(stream, " (%s)\n", str);
1676                     }
1677                 }
1678             }
1679         }
1680         /* It is possible to have a molecule without chains, in which case
1681          * residues in the molecule can be iterated through without going
1682          * through chains. */
1683         else
1684         {
1685             tng_molecule_num_residues_get(input, molecule, &nResidues);
1686             if (nResidues > 0)
1687             {
1688                 for (int64_t k = 0; k < nResidues; k++)
1689                 {
1690                     tng_molecule_residue_of_index_get(input, molecule, k, &residue);
1691                     tng_residue_name_get(input, residue, str, 256);
1692                     fprintf(stream, "\t\tResidue: %s\n", str);
1693                     tng_residue_num_atoms_get(input, residue, &nAtoms);
1694                     for (int64_t l = 0; l < nAtoms; l++)
1695                     {
1696                         tng_residue_atom_of_index_get(input, residue, l, &atom);
1697                         tng_atom_name_get(input, atom, str, 256);
1698                         fprintf(stream, "\t\t\tAtom: %s", str);
1699                         tng_atom_type_get(input, atom, str, 256);
1700                         fprintf(stream, " (%s)\n", str);
1701                     }
1702                 }
1703             }
1704             else
1705             {
1706                 tng_molecule_num_atoms_get(input, molecule, &nAtoms);
1707                 for (int64_t l = 0; l < nAtoms; l++)
1708                 {
1709                     tng_molecule_atom_of_index_get(input, molecule, l, &atom);
1710                     tng_atom_name_get(input, atom, str, 256);
1711                     fprintf(stream, "\t\t\tAtom: %s", str);
1712                     tng_atom_type_get(input, atom, str, 256);
1713                     fprintf(stream, " (%s)\n", str);
1714                 }
1715             }
1716         }
1717     }
1718
1719     tng_num_particles_get(input, &nAtoms);
1720     stat = tng_particle_data_vector_get(input, TNG_TRAJ_PARTIAL_CHARGES, &data, &nFramesRead,
1721                                         &strideLength, &nParticlesRead,
1722                                         &nValuesPerFrameRead, &datatype);
1723     if (stat == TNG_SUCCESS)
1724     {
1725         atomCharges.resize(nAtoms);
1726         convert_array_to_real_array(data,
1727                                     atomCharges.data(),
1728                                     1,
1729                                     nAtoms,
1730                                     1,
1731                                     datatype);
1732
1733         fprintf(stream, "Atom Charges (%d):\n", int(nAtoms));
1734         for (int64_t i = 0; i < nAtoms; i += 10)
1735         {
1736             fprintf(stream, "Atom Charges [%8d-]=[", int(i));
1737             for (int64_t j = 0; (j < 10 && i + j < nAtoms); j++)
1738             {
1739                 fprintf(stream, " %12.5e", atomCharges[i + j]);
1740             }
1741             fprintf(stream, "]\n");
1742         }
1743     }
1744
1745     stat = tng_particle_data_vector_get(input, TNG_TRAJ_MASSES, &data, &nFramesRead,
1746                                         &strideLength, &nParticlesRead,
1747                                         &nValuesPerFrameRead, &datatype);
1748     if (stat == TNG_SUCCESS)
1749     {
1750         atomMasses.resize(nAtoms);
1751         convert_array_to_real_array(data,
1752                                     atomMasses.data(),
1753                                     1,
1754                                     nAtoms,
1755                                     1,
1756                                     datatype);
1757
1758         fprintf(stream, "Atom Masses (%d):\n", int(nAtoms));
1759         for (int64_t i = 0; i < nAtoms; i += 10)
1760         {
1761             fprintf(stream, "Atom Masses [%8d-]=[", int(i));
1762             for (int64_t j = 0; (j < 10 && i + j < nAtoms); j++)
1763             {
1764                 fprintf(stream, " %12.5e", atomMasses[i + j]);
1765             }
1766             fprintf(stream, "]\n");
1767         }
1768     }
1769
1770     sfree(data);
1771 #else
1772     GMX_UNUSED_VALUE(gmx_tng_input);
1773     GMX_UNUSED_VALUE(stream);
1774 #endif
1775 }
1776
1777 gmx_bool gmx_get_tng_data_block_types_of_next_frame(gmx_tng_trajectory_t gmx_tng_input,
1778                                                     int                  frame,
1779                                                     int                  nRequestedIds,
1780                                                     int64_t             *requestedIds,
1781                                                     int64_t             *nextFrame,
1782                                                     int64_t             *nBlocks,
1783                                                     int64_t            **blockIds)
1784 {
1785 #if GMX_USE_TNG
1786     tng_function_status stat;
1787     tng_trajectory_t    input = gmx_tng_input->tng;
1788
1789     stat = tng_util_trajectory_next_frame_present_data_blocks_find(input, frame,
1790                                                                    nRequestedIds, requestedIds,
1791                                                                    nextFrame,
1792                                                                    nBlocks, blockIds);
1793
1794     if (stat == TNG_CRITICAL)
1795     {
1796         gmx_file("Cannot read TNG file. Cannot find data blocks of next frame.");
1797     }
1798     else if (stat == TNG_FAILURE)
1799     {
1800         return FALSE;
1801     }
1802     return TRUE;
1803 #else
1804     GMX_UNUSED_VALUE(gmx_tng_input);
1805     GMX_UNUSED_VALUE(frame);
1806     GMX_UNUSED_VALUE(nRequestedIds);
1807     GMX_UNUSED_VALUE(requestedIds);
1808     GMX_UNUSED_VALUE(nextFrame);
1809     GMX_UNUSED_VALUE(nBlocks);
1810     GMX_UNUSED_VALUE(blockIds);
1811     return FALSE;
1812 #endif
1813 }
1814
1815 gmx_bool gmx_get_tng_data_next_frame_of_block_type(gmx_tng_trajectory_t gmx_tng_input,
1816                                                    int64_t              blockId,
1817                                                    real               **values,
1818                                                    int64_t             *frameNumber,
1819                                                    double              *frameTime,
1820                                                    int64_t             *nValuesPerFrame,
1821                                                    int64_t             *nAtoms,
1822                                                    real                *prec,
1823                                                    char                *name,
1824                                                    int                  maxLen,
1825                                                    gmx_bool            *bOK)
1826 {
1827 #if GMX_USE_TNG
1828     tng_function_status stat;
1829     char                datatype = -1;
1830     int64_t             codecId;
1831     int                 blockDependency;
1832     void               *data = nullptr;
1833     double              localPrec;
1834     tng_trajectory_t    input = gmx_tng_input->tng;
1835
1836     stat = tng_data_block_name_get(input, blockId, name, maxLen);
1837     if (stat != TNG_SUCCESS)
1838     {
1839         gmx_file("Cannot read next frame of TNG file");
1840     }
1841     stat = tng_data_block_dependency_get(input, blockId, &blockDependency);
1842     if (stat != TNG_SUCCESS)
1843     {
1844         gmx_file("Cannot read next frame of TNG file");
1845     }
1846     if (blockDependency & TNG_PARTICLE_DEPENDENT)
1847     {
1848         tng_num_particles_get(input, nAtoms);
1849         stat = tng_util_particle_data_next_frame_read(input,
1850                                                       blockId,
1851                                                       &data,
1852                                                       &datatype,
1853                                                       frameNumber,
1854                                                       frameTime);
1855     }
1856     else
1857     {
1858         *nAtoms = 1; /* There are not actually any atoms, but it is used for
1859                         allocating memory */
1860         stat    = tng_util_non_particle_data_next_frame_read(input,
1861                                                              blockId,
1862                                                              &data,
1863                                                              &datatype,
1864                                                              frameNumber,
1865                                                              frameTime);
1866     }
1867     if (stat == TNG_CRITICAL)
1868     {
1869         gmx_file("Cannot read next frame of TNG file");
1870     }
1871     if (stat == TNG_FAILURE)
1872     {
1873         *bOK = TRUE;
1874         return FALSE;
1875     }
1876
1877     stat = tng_data_block_num_values_per_frame_get(input, blockId, nValuesPerFrame);
1878     if (stat != TNG_SUCCESS)
1879     {
1880         gmx_file("Cannot read next frame of TNG file");
1881     }
1882     srenew(*values, sizeof(real) * *nValuesPerFrame * *nAtoms);
1883     convert_array_to_real_array(data,
1884                                 *values,
1885                                 getDistanceScaleFactor(gmx_tng_input),
1886                                 *nAtoms,
1887                                 *nValuesPerFrame,
1888                                 datatype);
1889
1890     tng_util_frame_current_compression_get(input, blockId, &codecId, &localPrec);
1891
1892     /* This must be updated if/when more lossy compression methods are added */
1893     if (codecId != TNG_TNG_COMPRESSION)
1894     {
1895         *prec = -1.0;
1896     }
1897     else
1898     {
1899         *prec = localPrec;
1900     }
1901
1902     sfree(data);
1903
1904     *bOK = TRUE;
1905     return TRUE;
1906 #else
1907     GMX_UNUSED_VALUE(gmx_tng_input);
1908     GMX_UNUSED_VALUE(blockId);
1909     GMX_UNUSED_VALUE(values);
1910     GMX_UNUSED_VALUE(frameNumber);
1911     GMX_UNUSED_VALUE(frameTime);
1912     GMX_UNUSED_VALUE(nValuesPerFrame);
1913     GMX_UNUSED_VALUE(nAtoms);
1914     GMX_UNUSED_VALUE(prec);
1915     GMX_UNUSED_VALUE(name);
1916     GMX_UNUSED_VALUE(maxLen);
1917     GMX_UNUSED_VALUE(bOK);
1918     return FALSE;
1919 #endif
1920 }
1921
1922 int gmx_tng_get_box_output_interval(gmx_tng_trajectory_t gmx_tng)
1923 {
1924 #if GMX_USE_TNG
1925     return gmx_tng->boxOutputInterval;
1926 #else
1927     GMX_UNUSED_VALUE(gmx_tng);
1928 #endif
1929 }
1930
1931 int gmx_tng_get_lambda_output_interval(gmx_tng_trajectory_t gmx_tng)
1932 {
1933 #if GMX_USE_TNG
1934     return gmx_tng->lambdaOutputInterval;
1935 #else
1936     GMX_UNUSED_VALUE(gmx_tng);
1937 #endif
1938 }