src/gromacs/fileio/tngio.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include "gmxpre.h"
  36
  37 #include "tngio.h"
  38
  39 #include "config.h"
  40
  41 #ifdef HAVE_UNISTD_H
  42 #include <unistd.h>
  43 #endif
  44
  45 #ifdef GMX_USE_TNG
  46 #include "tng/tng_io.h"
  47 #endif
  48
  49 #include "gromacs/legacyheaders/copyrite.h"
  50 #include "gromacs/legacyheaders/types/ifunc.h"
  51
  52 #include "gromacs/fileio/gmxfio.h"
  53 #include "gromacs/math/units.h"
  54 #include "gromacs/math/utilities.h"
  55 #include "gromacs/topology/topology.h"
  56 #include "gromacs/utility/basenetwork.h"
  57 #include "gromacs/utility/common.h"
  58 #include "gromacs/utility/fatalerror.h"
  59 #include "gromacs/utility/gmxassert.h"
  60 #include "gromacs/utility/programcontext.h"
  61
  62 static const char *modeToVerb(char mode)
  63 {
  64     const char *p;
  65     switch (mode)
  66     {
  67         case 'r':
  68             p = "reading";
  69             break;
  70         case 'w':
  71             p = "writing";
  72             break;
  73         case 'a':
  74             p = "appending";
  75             break;
  76         default:
  77             gmx_fatal(FARGS, "Invalid file opening mode %c", mode);
  78             p = "";
  79             break;
  80     }
  81     return p;
  82 }
  83
  84 void gmx_tng_open(const char       *filename,
  85                   char              mode,
  86                   tng_trajectory_t *tng)
  87 {
  88 #ifdef GMX_USE_TNG
  89     /* First check whether we have to make a backup,
  90      * only for writing, not for read or append.
  91      */
  92     if (mode == 'w')
  93     {
  94 #ifndef GMX_FAHCORE
  95         /* only make backups for normal gromacs */
  96         make_backup(filename);
  97 #endif
  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                   "%s while opening %s for %s",
 110                   gmx_strerror("file"),
 111                   filename,
 112                   modeToVerb(mode));
 113     }
 114
 115     if (mode == 'w' || mode == 'a')
 116     {
 117         /* FIXME in TNG: When adding data to the header, subsequent blocks might get
 118          * overwritten. This could be solved by moving the first trajectory
 119          * frame set(s) to the end of the file. Could that cause other problems,
 120          * e.g. when continuing a simulation? */
 121         char hostname[256];
 122         gmx_gethostname(hostname, 256);
 123         if (mode == 'w')
 124         {
 125             tng_first_computer_name_set(*tng, hostname);
 126         }
 127 /* TODO: This should be implemented when the above fixme is done (adding data to
 128  * the header). */
 129 //         else
 130 //         {
 131 //             tng_last_computer_name_set(*tng, hostname);
 132 //         }
 133
 134         char        programInfo[256];
 135         const char *precisionString = "";
 136 #ifdef GMX_DOUBLE
 137         precisionString = " (double precision)";
 138 #endif
 139         sprintf(programInfo, "%.100s, %.128s%.24s",
 140                 gmx::getProgramContext().displayName(),
 141                 GromacsVersion(), precisionString);
 142         if (mode == 'w')
 143         {
 144             tng_first_program_name_set(*tng, programInfo);
 145         }
 146 /* TODO: This should be implemented when the above fixme is done (adding data to
 147  * the header). */
 148 //         else
 149 //         {
 150 //             tng_last_program_name_set(*tng, programInfo);
 151 //         }
 152
 153 #if defined(HAVE_UNISTD_H) && !defined(__MINGW32__)
 154         char username[256];
 155         if (!getlogin_r(username, 256))
 156         {
 157             if (mode == 'w')
 158             {
 159                 tng_first_user_name_set(*tng, username);
 160             }
 161         }
 162 /* TODO: This should be implemented when the above fixme is done (adding data to
 163  * the header). */
 164 //         else
 165 //         {
 166 //             tng_last_user_name_set(*tng, username);
 167 //         }
 168 #endif
 169     }
 170 #else
 171     gmx_file("GROMACS was compiled without TNG support, cannot handle this file type");
 172     GMX_UNUSED_VALUE(filename);
 173     GMX_UNUSED_VALUE(mode);
 174     GMX_UNUSED_VALUE(tng);
 175 #endif
 176 }
 177
 178 void gmx_tng_close(tng_trajectory_t *tng)
 179 {
 180     /* We have to check that tng is set because
 181      * tng_util_trajectory_close wants to return a NULL in it, and
 182      * gives a fatal error if it is NULL. */
 183 #ifdef GMX_USE_TNG
 184     if (tng)
 185     {
 186         tng_util_trajectory_close(tng);
 187     }
 188 #else
 189     GMX_UNUSED_VALUE(tng);
 190 #endif
 191 }
 192
 193 #ifdef GMX_USE_TNG
 194 static void addTngMoleculeFromTopology(tng_trajectory_t     tng,
 195                                        const char          *moleculeName,
 196                                        const t_atoms       *atoms,
 197                                        gmx_int64_t          numMolecules,
 198                                        tng_molecule_t      *tngMol)
 199 {
 200     tng_chain_t      tngChain = NULL;
 201     tng_residue_t    tngRes   = NULL;
 202
 203     if (tng_molecule_add(tng, moleculeName, tngMol) != TNG_SUCCESS)
 204     {
 205         gmx_file("Cannot add molecule to TNG molecular system.");
 206     }
 207
 208     /* FIXME: The TNG atoms should contain mass and atomB info (for free
 209      * energy calculations), i.e. in when it's available in TNG (2.0). */
 210     for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++)
 211     {
 212         const t_atom *at = &atoms->atom[atomIndex];
 213         /* FIXME: Currently the TNG API can only add atoms belonging to a
 214          * residue and chain. Wait for TNG 2.0*/
 215         if (atoms->nres > 0)
 216         {
 217             const t_resinfo *resInfo        = &atoms->resinfo[at->resind];
 218             char             chainName[2]   = {resInfo->chainid, 0};
 219             tng_atom_t       tngAtom        = NULL;
 220             t_atom          *prevAtom;
 221
 222             if (atomIndex > 0)
 223             {
 224                 prevAtom = &atoms->atom[atomIndex - 1];
 225             }
 226             else
 227             {
 228                 prevAtom = 0;
 229             }
 230
 231             /* If this is the first atom or if the residue changed add the
 232              * residue to the TNG molecular system. */
 233             if (!prevAtom || resInfo != &atoms->resinfo[prevAtom->resind])
 234             {
 235                 /* If this is the first atom or if the chain changed add
 236                  * the chain to the TNG molecular system. */
 237                 if (!prevAtom || resInfo->chainid !=
 238                     atoms->resinfo[prevAtom->resind].chainid)
 239                 {
 240                     tng_molecule_chain_add(tng, *tngMol, chainName,
 241                                            &tngChain);
 242                 }
 243                 /* FIXME: When TNG supports both residue index and residue
 244                  * number the latter should be used. Wait for TNG 2.0*/
 245                 tng_chain_residue_add(tng, tngChain, *resInfo->name, &tngRes);
 246             }
 247             tng_residue_atom_add(tng, tngRes, *(atoms->atomname[atomIndex]), *(atoms->atomtype[atomIndex]), &tngAtom);
 248         }
 249     }
 250     tng_molecule_cnt_set(tng, *tngMol, numMolecules);
 251 }
 252
 253 void gmx_tng_add_mtop(tng_trajectory_t  tng,
 254                       const gmx_mtop_t *mtop)
 255 {
 256     int                  i, j;
 257     const t_ilist       *ilist;
 258     tng_bond_t           tngBond;
 259
 260     if (!mtop)
 261     {
 262         /* No topology information available to add. */
 263         return;
 264     }
 265
 266     for (int molIndex = 0; molIndex < mtop->nmolblock; molIndex++)
 267     {
 268         tng_molecule_t       tngMol  = NULL;
 269         const gmx_moltype_t *molType =
 270             &mtop->moltype[mtop->molblock[molIndex].type];
 271
 272         /* Add a molecule to the TNG trajectory with the same name as the
 273          * current molecule. */
 274         addTngMoleculeFromTopology(tng,
 275                                    *(molType->name),
 276                                    &molType->atoms,
 277                                    mtop->molblock[molIndex].nmol,
 278                                    &tngMol);
 279
 280         /* Bonds have to be deduced from interactions (constraints etc). Different
 281          * interactions have different sets of parameters. */
 282         /* Constraints are specified using two atoms */
 283         for (i = 0; i < F_NRE; i++)
 284         {
 285             if (IS_CHEMBOND(i))
 286             {
 287                 ilist = &molType->ilist[i];
 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                         j += 3;
 295                     }
 296                 }
 297             }
 298         }
 299         /* Settle is described using three atoms */
 300         ilist = &molType->ilist[F_SETTLE];
 301         if (ilist)
 302         {
 303             j = 1;
 304             while (j < ilist->nr)
 305             {
 306                 tng_molecule_bond_add(tng, tngMol, ilist->iatoms[j], ilist->iatoms[j+1], &tngBond);
 307                 tng_molecule_bond_add(tng, tngMol, ilist->iatoms[j], ilist->iatoms[j+2], &tngBond);
 308                 j += 4;
 309             }
 310         }
 311     }
 312 }
 313
 314 /*! \libinternal \brief Compute greatest common divisor of n1 and n2
 315  * if they are positive.
 316  *
 317  * If only one of n1 and n2 is positive, then return it.
 318  * If neither n1 or n2 is positive, then return -1. */
 319 static int
 320 greatest_common_divisor_if_positive(int n1, int n2)
 321 {
 322     if (0 >= n1)
 323     {
 324         return (0 >= n2) ? -1 : n2;
 325     }
 326     if (0 >= n2)
 327     {
 328         return n1;
 329     }
 330
 331     /* We have a non-trivial greatest common divisor to compute. */
 332     return gmx_greatest_common_divisor(n1, n2);
 333 }
 334
 335 /* By default try to write 100 frames (of actual output) in each frame set.
 336  * This number is the number of outputs of the most frequently written data
 337  * type per frame set.
 338  * TODO for 5.1: Verify that 100 frames per frame set is efficient for most
 339  * setups regarding compression efficiency and compression time. Make this
 340  * a hidden command-line option? */
 341 const int defaultFramesPerFrameSet = 100;
 342
 343 /*! \libinternal \brief  Set the number of frames per frame
 344  * set according to output intervals.
 345  * The default is that 100 frames are written of the data
 346  * that is written most often. */
 347 static void tng_set_frames_per_frame_set(tng_trajectory_t  tng,
 348                                          const gmx_bool    bUseLossyCompression,
 349                                          const t_inputrec *ir)
 350 {
 351     int     gcd = -1;
 352
 353     /* Set the number of frames per frame set to contain at least
 354      * defaultFramesPerFrameSet of the lowest common denominator of
 355      * the writing interval of positions and velocities. */
 356     /* FIXME after 5.0: consider nstenergy also? */
 357     if (bUseLossyCompression)
 358     {
 359         gcd = ir->nstxout_compressed;
 360     }
 361     else
 362     {
 363         gcd = greatest_common_divisor_if_positive(ir->nstxout, ir->nstvout);
 364         gcd = greatest_common_divisor_if_positive(gcd, ir->nstfout);
 365     }
 366     if (0 >= gcd)
 367     {
 368         return;
 369     }
 370
 371     tng_num_frames_per_frame_set_set(tng, gcd * defaultFramesPerFrameSet);
 372 }
 373
 374 /*! \libinternal \brief Set the data-writing intervals, and number of
 375  * frames per frame set */
 376 static void set_writing_intervals(tng_trajectory_t  tng,
 377                                   const gmx_bool    bUseLossyCompression,
 378                                   const t_inputrec *ir)
 379 {
 380     /* Define pointers to specific writing functions depending on if we
 381      * write float or double data */
 382     typedef tng_function_status (*set_writing_interval_func_pointer)(tng_trajectory_t,
 383                                                                      const gmx_int64_t,
 384                                                                      const gmx_int64_t,
 385                                                                      const gmx_int64_t,
 386                                                                      const char*,
 387                                                                      const char,
 388                                                                      const char);
 389 #ifdef GMX_DOUBLE
 390     set_writing_interval_func_pointer set_writing_interval = tng_util_generic_write_interval_double_set;
 391 #else
 392     set_writing_interval_func_pointer set_writing_interval = tng_util_generic_write_interval_set;
 393 #endif
 394     int  xout, vout, fout;
 395 //     int  gcd = -1, lowest = -1;
 396     char compression;
 397
 398     tng_set_frames_per_frame_set(tng, bUseLossyCompression, ir);
 399
 400     if (bUseLossyCompression)
 401     {
 402         xout        = ir->nstxout_compressed;
 403         vout        = 0;
 404         fout        = 0;
 405         compression = TNG_TNG_COMPRESSION;
 406     }
 407     else
 408     {
 409         xout        = ir->nstxout;
 410         vout        = ir->nstvout;
 411         fout        = ir->nstfout;
 412         compression = TNG_GZIP_COMPRESSION;
 413     }
 414     if (xout)
 415     {
 416         set_writing_interval(tng, xout, 3, TNG_TRAJ_POSITIONS,
 417                              "POSITIONS", TNG_PARTICLE_BLOCK_DATA,
 418                              compression);
 419         /* The design of TNG makes it awkward to try to write a box
 420          * with multiple periodicities, which might be co-prime. Since
 421          * the use cases for the box with a frame consisting only of
 422          * velocities seem low, for now we associate box writing with
 423          * position writing. */
 424         set_writing_interval(tng, xout, 9, TNG_TRAJ_BOX_SHAPE,
 425                              "BOX SHAPE", TNG_NON_PARTICLE_BLOCK_DATA,
 426                              TNG_GZIP_COMPRESSION);
 427         /* TODO: if/when we write energies to TNG also, reconsider how
 428          * and when box information is written, because GROMACS
 429          * behaviour pre-5.0 was to write the box with every
 430          * trajectory frame and every energy frame, and probably
 431          * people depend on this. */
 432
 433         /* TODO: If we need to write lambda values at steps when
 434          * positions (or other data) are not also being written, then
 435          * code in mdoutf.c will need to match however that is
 436          * organized here. */
 437         set_writing_interval(tng, xout, 1, TNG_GMX_LAMBDA,
 438                              "LAMBDAS", TNG_NON_PARTICLE_BLOCK_DATA,
 439                              TNG_GZIP_COMPRESSION);
 440
 441         /* FIXME: gcd and lowest currently not used. */
 442 //         gcd = greatest_common_divisor_if_positive(gcd, xout);
 443 //         if (lowest < 0 || xout < lowest)
 444 //         {
 445 //             lowest = xout;
 446 //         }
 447     }
 448     if (vout)
 449     {
 450         set_writing_interval(tng, ir->nstvout, 3, TNG_TRAJ_VELOCITIES,
 451                              "VELOCITIES", TNG_PARTICLE_BLOCK_DATA,
 452                              compression);
 453
 454         /* FIXME: gcd and lowest currently not used. */
 455 //         gcd = greatest_common_divisor_if_positive(gcd, vout);
 456 //         if (lowest < 0 || vout < lowest)
 457 //         {
 458 //             lowest = vout;
 459 //         }
 460     }
 461     if (fout)
 462     {
 463         set_writing_interval(tng, ir->nstfout, 3, TNG_TRAJ_FORCES,
 464                              "FORCES", TNG_PARTICLE_BLOCK_DATA,
 465                              TNG_GZIP_COMPRESSION);
 466
 467         /* FIXME: gcd and lowest currently not used. */
 468 //         gcd = greatest_common_divisor_if_positive(gcd, fout);
 469 //         if (lowest < 0 || fout < lowest)
 470 //         {
 471 //             lowest = fout;
 472 //         }
 473     }
 474     /* FIXME: See above. gcd interval for lambdas is disabled. */
 475 //     if (gcd > 0)
 476 //     {
 477 //         /* Lambdas written at an interval of the lowest common denominator
 478 //          * of other output */
 479 //         set_writing_interval(tng, gcd, 1, TNG_GMX_LAMBDA,
 480 //                                  "LAMBDAS", TNG_NON_PARTICLE_BLOCK_DATA,
 481 //                                  TNG_GZIP_COMPRESSION);
 482 //
 483 //         if (gcd < lowest / 10)
 484 //         {
 485 //             gmx_warning("The lowest common denominator of trajectory output is "
 486 //                         "every %d step(s), whereas the shortest output interval "
 487 //                         "is every %d steps.", gcd, lowest);
 488 //         }
 489 //     }
 490 }
 491 #endif
 492
 493 void gmx_tng_prepare_md_writing(tng_trajectory_t  tng,
 494                                 const gmx_mtop_t *mtop,
 495                                 const t_inputrec *ir)
 496 {
 497 #ifdef GMX_USE_TNG
 498     gmx_tng_add_mtop(tng, mtop);
 499     set_writing_intervals(tng, FALSE, ir);
 500     tng_time_per_frame_set(tng, ir->delta_t * PICO);
 501 #else
 502     GMX_UNUSED_VALUE(tng);
 503     GMX_UNUSED_VALUE(mtop);
 504     GMX_UNUSED_VALUE(ir);
 505 #endif
 506 }
 507
 508 #ifdef GMX_USE_TNG
 509 /* Check if all atoms in the molecule system are selected
 510  * by a selection group of type specified by gtype. */
 511 static gmx_bool all_atoms_selected(const gmx_mtop_t *mtop,
 512                                    const int         gtype)
 513 {
 514     const gmx_moltype_t     *molType;
 515     const t_atoms           *atoms;
 516
 517     /* Iterate through all atoms in the molecule system and
 518      * check if they belong to a selection group of the
 519      * requested type. */
 520     for (int molIndex = 0, i = 0; molIndex < mtop->nmoltype; molIndex++)
 521     {
 522         molType = &mtop->moltype[mtop->molblock[molIndex].type];
 523
 524         atoms = &molType->atoms;
 525
 526         for (int j = 0; j < mtop->molblock[molIndex].nmol; j++)
 527         {
 528             for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++, i++)
 529             {
 530                 if (ggrpnr(&mtop->groups, gtype, i) != 0)
 531                 {
 532                     return FALSE;
 533                 }
 534             }
 535         }
 536     }
 537     return TRUE;
 538 }
 539
 540 /* Create TNG molecules which will represent each of the selection
 541  * groups for writing. But do that only if there is actually a
 542  * specified selection group and it is not the whole system.
 543  * TODO: Currently the only selection that is taken into account
 544  * is egcCompressedX, but other selections should be added when
 545  * e.g. writing energies is implemented.
 546  */
 547 static void add_selection_groups(tng_trajectory_t  tng,
 548                                  const gmx_mtop_t *mtop)
 549 {
 550     const gmx_moltype_t     *molType;
 551     const t_atoms           *atoms;
 552     const t_atom            *at;
 553     const t_resinfo         *resInfo;
 554     const t_ilist           *ilist;
 555     int                      nameIndex;
 556     int                      atom_offset = 0;
 557     tng_molecule_t           mol, iterMol;
 558     tng_chain_t              chain;
 559     tng_residue_t            res;
 560     tng_atom_t               atom;
 561     tng_bond_t               tngBond;
 562     gmx_int64_t              nMols;
 563     char                    *groupName;
 564
 565     /* TODO: When the TNG molecules block is more flexible TNG selection
 566      * groups should not need all atoms specified. It should be possible
 567      * just to specify what molecules are selected (and/or which atoms in
 568      * the molecule) and how many (if applicable). */
 569
 570     /* If no atoms are selected we do not need to create a
 571      * TNG selection group molecule. */
 572     if (mtop->groups.ngrpnr[egcCompressedX] == 0)
 573     {
 574         return;
 575     }
 576
 577     /* If all atoms are selected we do not have to create a selection
 578      * group molecule in the TNG molecule system. */
 579     if (all_atoms_selected(mtop, egcCompressedX))
 580     {
 581         return;
 582     }
 583
 584     /* The name of the TNG molecule containing the selection group is the
 585      * same as the name of the selection group. */
 586     nameIndex = *mtop->groups.grps[egcCompressedX].nm_ind;
 587     groupName = *mtop->groups.grpname[nameIndex];
 588
 589     tng_molecule_alloc(tng, &mol);
 590     tng_molecule_name_set(tng, mol, groupName);
 591     tng_molecule_chain_add(tng, mol, "", &chain);
 592     for (int molIndex = 0, i = 0; molIndex < mtop->nmoltype; molIndex++)
 593     {
 594         molType = &mtop->moltype[mtop->molblock[molIndex].type];
 595
 596         atoms = &molType->atoms;
 597
 598         for (int j = 0; j < mtop->molblock[molIndex].nmol; j++)
 599         {
 600             bool bAtomsAdded = FALSE;
 601             for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++, i++)
 602             {
 603                 char *res_name;
 604                 int   res_id;
 605
 606                 if (ggrpnr(&mtop->groups, egcCompressedX, i) != 0)
 607                 {
 608                     continue;
 609                 }
 610                 at = &atoms->atom[atomIndex];
 611                 if (atoms->nres > 0)
 612                 {
 613                     resInfo = &atoms->resinfo[at->resind];
 614                     /* FIXME: When TNG supports both residue index and residue
 615                      * number the latter should be used. */
 616                     res_name = *resInfo->name;
 617                     res_id   = at->resind + 1;
 618                 }
 619                 else
 620                 {
 621                     res_name = (char *)"";
 622                     res_id   = 0;
 623                 }
 624                 if (tng_chain_residue_find(tng, chain, res_name, res_id, &res)
 625                     != TNG_SUCCESS)
 626                 {
 627                     /* Since there is ONE chain for selection groups do not keep the
 628                      * original residue IDs - otherwise there might be conflicts. */
 629                     tng_chain_residue_add(tng, chain, res_name, &res);
 630                 }
 631                 tng_residue_atom_w_id_add(tng, res, *(atoms->atomname[atomIndex]),
 632                                           *(atoms->atomtype[atomIndex]),
 633                                           atom_offset + atomIndex, &atom);
 634                 bAtomsAdded = TRUE;
 635             }
 636             /* Add bonds. */
 637             if (bAtomsAdded)
 638             {
 639                 for (int k = 0; k < F_NRE; k++)
 640                 {
 641                     if (IS_CHEMBOND(k))
 642                     {
 643                         ilist = &molType->ilist[k];
 644                         if (ilist)
 645                         {
 646                             int l = 1;
 647                             while (l < ilist->nr)
 648                             {
 649                                 int atom1, atom2;
 650                                 atom1 = ilist->iatoms[l] + atom_offset;
 651                                 atom2 = ilist->iatoms[l+1] + atom_offset;
 652                                 if (ggrpnr(&mtop->groups, egcCompressedX, atom1) == 0 &&
 653                                     ggrpnr(&mtop->groups, egcCompressedX, atom2) == 0)
 654                                 {
 655                                     tng_molecule_bond_add(tng, mol, ilist->iatoms[l],
 656                                                           ilist->iatoms[l+1], &tngBond);
 657                                 }
 658                                 l += 3;
 659                             }
 660                         }
 661                     }
 662                 }
 663                 /* Settle is described using three atoms */
 664                 ilist = &molType->ilist[F_SETTLE];
 665                 if (ilist)
 666                 {
 667                     int l = 1;
 668                     while (l < ilist->nr)
 669                     {
 670                         int atom1, atom2, atom3;
 671                         atom1 = ilist->iatoms[l] + atom_offset;
 672                         atom2 = ilist->iatoms[l+1] + atom_offset;
 673                         atom3 = ilist->iatoms[l+2] + atom_offset;
 674                         if (ggrpnr(&mtop->groups, egcCompressedX, atom1) == 0)
 675                         {
 676                             if (ggrpnr(&mtop->groups, egcCompressedX, atom2) == 0)
 677                             {
 678                                 tng_molecule_bond_add(tng, mol, atom1,
 679                                                       atom2, &tngBond);
 680                             }
 681                             if (ggrpnr(&mtop->groups, egcCompressedX, atom3) == 0)
 682                             {
 683                                 tng_molecule_bond_add(tng, mol, atom1,
 684                                                       atom3, &tngBond);
 685                             }
 686                         }
 687                         l += 4;
 688                     }
 689                 }
 690             }
 691             atom_offset += atoms->nr;
 692         }
 693     }
 694     tng_molecule_existing_add(tng, &mol);
 695     tng_molecule_cnt_set(tng, mol, 1);
 696     tng_num_molecule_types_get(tng, &nMols);
 697     for (gmx_int64_t k = 0; k < nMols; k++)
 698     {
 699         tng_molecule_of_index_get(tng, k, &iterMol);
 700         if (iterMol == mol)
 701         {
 702             continue;
 703         }
 704         tng_molecule_cnt_set(tng, iterMol, 0);
 705     }
 706 }
 707 #endif
 708
 709 void gmx_tng_set_compression_precision(tng_trajectory_t tng,
 710                                        real             prec)
 711 {
 712 #ifdef GMX_USE_TNG
 713     tng_compression_precision_set(tng, prec);
 714 #else
 715     GMX_UNUSED_VALUE(tng);
 716     GMX_UNUSED_VALUE(prec);
 717 #endif
 718 }
 719
 720 void gmx_tng_prepare_low_prec_writing(tng_trajectory_t  tng,
 721                                       const gmx_mtop_t *mtop,
 722                                       const t_inputrec *ir)
 723 {
 724 #ifdef GMX_USE_TNG
 725     gmx_tng_add_mtop(tng, mtop);
 726     add_selection_groups(tng, mtop);
 727     set_writing_intervals(tng, TRUE, ir);
 728     tng_time_per_frame_set(tng, ir->delta_t * PICO);
 729     gmx_tng_set_compression_precision(tng, ir->x_compression_precision);
 730 #else
 731     GMX_UNUSED_VALUE(tng);
 732     GMX_UNUSED_VALUE(mtop);
 733     GMX_UNUSED_VALUE(ir);
 734 #endif
 735 }
 736
 737 void gmx_fwrite_tng(tng_trajectory_t tng,
 738                     const gmx_bool   bUseLossyCompression,
 739                     int              step,
 740                     real             elapsedPicoSeconds,
 741                     real             lambda,
 742                     const rvec      *box,
 743                     int              nAtoms,
 744                     const rvec      *x,
 745                     const rvec      *v,
 746                     const rvec      *f)
 747 {
 748 #ifdef GMX_USE_TNG
 749     typedef tng_function_status (*write_data_func_pointer)(tng_trajectory_t,
 750                                                            const gmx_int64_t,
 751                                                            const double,
 752                                                            const real*,
 753                                                            const gmx_int64_t,
 754                                                            const gmx_int64_t,
 755                                                            const char*,
 756                                                            const char,
 757                                                            const char);
 758 #ifdef GMX_DOUBLE
 759     static write_data_func_pointer           write_data           = tng_util_generic_with_time_double_write;
 760 #else
 761     static write_data_func_pointer           write_data           = tng_util_generic_with_time_write;
 762 #endif
 763     double                                   elapsedSeconds = elapsedPicoSeconds * PICO;
 764     gmx_int64_t                              nParticles;
 765     char                                     compression;
 766
 767
 768     if (!tng)
 769     {
 770         /* This function might get called when the type of the
 771            compressed trajectory is actually XTC. So we exit and move
 772            on. */
 773         return;
 774     }
 775
 776     tng_num_particles_get(tng, &nParticles);
 777     if (nAtoms != (int)nParticles)
 778     {
 779         tng_implicit_num_particles_set(tng, nAtoms);
 780     }
 781
 782     if (bUseLossyCompression)
 783     {
 784         compression = TNG_TNG_COMPRESSION;
 785     }
 786     else
 787     {
 788         compression = TNG_GZIP_COMPRESSION;
 789     }
 790
 791     /* The writing is done using write_data, which writes float or double
 792      * depending on the GROMACS compilation. */
 793     if (x)
 794     {
 795         GMX_ASSERT(box, "Need a non-NULL box if positions are written");
 796
 797         if (write_data(tng, step, elapsedSeconds,
 798                        reinterpret_cast<const real *>(x),
 799                        3, TNG_TRAJ_POSITIONS, "POSITIONS",
 800                        TNG_PARTICLE_BLOCK_DATA,
 801                        compression) != TNG_SUCCESS)
 802         {
 803             gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 804         }
 805         /* TNG-MF1 compression only compresses positions and velocities. Use lossless
 806          * compression for box shape regardless of output mode */
 807         if (write_data(tng, step, elapsedSeconds,
 808                        reinterpret_cast<const real *>(box),
 809                        9, TNG_TRAJ_BOX_SHAPE, "BOX SHAPE",
 810                        TNG_NON_PARTICLE_BLOCK_DATA,
 811                        TNG_GZIP_COMPRESSION) != TNG_SUCCESS)
 812         {
 813             gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 814         }
 815     }
 816
 817     if (v)
 818     {
 819         if (write_data(tng, step, elapsedSeconds,
 820                        reinterpret_cast<const real *>(v),
 821                        3, TNG_TRAJ_VELOCITIES, "VELOCITIES",
 822                        TNG_PARTICLE_BLOCK_DATA,
 823                        compression) != TNG_SUCCESS)
 824         {
 825             gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 826         }
 827     }
 828
 829     if (f)
 830     {
 831         /* TNG-MF1 compression only compresses positions and velocities. Use lossless
 832          * compression for forces regardless of output mode */
 833         if (write_data(tng, step, elapsedSeconds,
 834                        reinterpret_cast<const real *>(f),
 835                        3, TNG_TRAJ_FORCES, "FORCES",
 836                        TNG_PARTICLE_BLOCK_DATA,
 837                        TNG_GZIP_COMPRESSION) != TNG_SUCCESS)
 838         {
 839             gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 840         }
 841     }
 842
 843     /* TNG-MF1 compression only compresses positions and velocities. Use lossless
 844      * compression for lambdas regardless of output mode */
 845     if (write_data(tng, step, elapsedSeconds,
 846                    reinterpret_cast<const real *>(&lambda),
 847                    1, TNG_GMX_LAMBDA, "LAMBDAS",
 848                    TNG_NON_PARTICLE_BLOCK_DATA,
 849                    TNG_GZIP_COMPRESSION) != TNG_SUCCESS)
 850     {
 851         gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 852     }
 853 #else
 854     GMX_UNUSED_VALUE(tng);
 855     GMX_UNUSED_VALUE(bUseLossyCompression);
 856     GMX_UNUSED_VALUE(step);
 857     GMX_UNUSED_VALUE(elapsedPicoSeconds);
 858     GMX_UNUSED_VALUE(lambda);
 859     GMX_UNUSED_VALUE(box);
 860     GMX_UNUSED_VALUE(nAtoms);
 861     GMX_UNUSED_VALUE(x);
 862     GMX_UNUSED_VALUE(v);
 863     GMX_UNUSED_VALUE(f);
 864 #endif
 865 }
 866
 867 void fflush_tng(tng_trajectory_t tng)
 868 {
 869 #ifdef GMX_USE_TNG
 870     if (!tng)
 871     {
 872         return;
 873     }
 874     tng_frame_set_premature_write(tng, TNG_USE_HASH);
 875 #else
 876     GMX_UNUSED_VALUE(tng);
 877 #endif
 878 }
 879
 880 float gmx_tng_get_time_of_final_frame(tng_trajectory_t tng)
 881 {
 882 #ifdef GMX_USE_TNG
 883     gmx_int64_t nFrames;
 884     double      time;
 885     float       fTime;
 886
 887     tng_num_frames_get(tng, &nFrames);
 888     tng_util_time_of_frame_get(tng, nFrames - 1, &time);
 889
 890     fTime = time / PICO;
 891     return fTime;
 892 #else
 893     GMX_UNUSED_VALUE(tng);
 894     return -1.0;
 895 #endif
 896 }