src/gromacs/fileio/tngio.cpp

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