src/gromacs/fileio/tngio.cpp

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