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