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