src/gromacs/fileio/tngio.cpp

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