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