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