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