src/gromacs/fileio/tngio.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 #include "tngio.h"
  36
  37 #ifdef HAVE_CONFIG_H
  38 #include <config.h>
  39 #endif
  40
  41 #ifdef HAVE_UNISTD_H
  42 #include <unistd.h>
  43 #endif
  44
  45 #ifdef GMX_USE_TNG
  46 #include "../../external/tng_io/include/tng_io.h"
  47 #endif
  48
  49 #include "gromacs/legacyheaders/copyrite.h"
  50 #include "gromacs/legacyheaders/gmx_fatal.h"
  51 #include "gromacs/legacyheaders/main.h"
  52 #include "gromacs/legacyheaders/physics.h"
  53 #include "gromacs/math/utilities.h"
  54 #include "gromacs/utility/common.h"
  55 #include "gromacs/utility/gmxassert.h"
  56 #include "gromacs/utility/programcontext.h"
  57 #include "gmxfio.h"
  58
  59 static const char *modeToVerb(char mode)
  60 {
  61     switch (mode)
  62     {
  63         case 'r':
  64             return "reading";
  65             break;
  66         case 'w':
  67             return "writing";
  68             break;
  69         case 'a':
  70             return "appending";
  71             break;
  72         default:
  73             gmx_fatal(FARGS, "Invalid file opening mode %c", mode);
  74             return "";
  75     }
  76 }
  77
  78 void gmx_tng_open(const char       *filename,
  79                   char              mode,
  80                   tng_trajectory_t *tng)
  81 {
  82 #ifdef GMX_USE_TNG
  83     /* First check whether we have to make a backup,
  84      * only for writing, not for read or append.
  85      */
  86     if (mode == 'w')
  87     {
  88 #ifndef GMX_FAHCORE
  89         /* only make backups for normal gromacs */
  90         make_backup(filename);
  91 #endif
  92     }
  93
  94     /* tng must not be pointing at already allocated memory.
  95      * Memory will be allocated by tng_util_trajectory_open() and must
  96      * later on be freed by tng_util_trajectory_close(). */
  97     if (TNG_SUCCESS != tng_util_trajectory_open(filename, mode, tng))
  98     {
  99         /* TNG does return more than one degree of error, but there is
 100            no use case for GROMACS handling the non-fatal errors
 101            gracefully. */
 102         gmx_fatal(FARGS,
 103                   "%s while opening %s for %s",
 104                   gmx_strerror("file"),
 105                   filename,
 106                   modeToVerb(mode));
 107     }
 108
 109     if (mode == 'w' || mode == 'a')
 110     {
 111         /* FIXME in TNG: When adding data to the header, subsequent blocks might get
 112          * overwritten. This could be solved by moving the first trajectory
 113          * frame set(s) to the end of the file. Could that cause other problems,
 114          * e.g. when continuing a simulation? */
 115         char hostname[256];
 116         gmx_gethostname(hostname, 256);
 117         if (mode == 'w')
 118         {
 119             tng_first_computer_name_set(*tng, hostname);
 120         }
 121 /* TODO: This should be implemented when the above fixme is done (adding data to
 122  * the header). */
 123 //         else
 124 //         {
 125 //             tng_last_computer_name_set(*tng, hostname);
 126 //         }
 127
 128         char        programInfo[256];
 129         const char *precisionString = "";
 130 #ifdef GMX_DOUBLE
 131         precisionString = " (double precision)";
 132 #endif
 133         sprintf(programInfo, "%.100s, %.128s%.24s",
 134                 gmx::getProgramContext().displayName(),
 135                 GromacsVersion(), precisionString);
 136         if (mode == 'w')
 137         {
 138             tng_first_program_name_set(*tng, programInfo);
 139         }
 140 /* TODO: This should be implemented when the above fixme is done (adding data to
 141  * the header). */
 142 //         else
 143 //         {
 144 //             tng_last_program_name_set(*tng, programInfo);
 145 //         }
 146
 147 #ifdef HAVE_UNISTD_H
 148         char username[256];
 149         getlogin_r(username, 256);
 150         if (mode == 'w')
 151         {
 152             tng_first_user_name_set(*tng, username);
 153         }
 154 /* TODO: This should be implemented when the above fixme is done (adding data to
 155  * the header). */
 156 //         else
 157 //         {
 158 //             tng_last_user_name_set(*tng, username);
 159 //         }
 160 #endif
 161     }
 162 #else
 163     gmx_file("GROMACS was compiled without TNG support, cannot handle this file type");
 164     GMX_UNUSED_VALUE(filename);
 165     GMX_UNUSED_VALUE(mode);
 166     GMX_UNUSED_VALUE(tng);
 167 #endif
 168 }
 169
 170 void gmx_tng_close(tng_trajectory_t *tng)
 171 {
 172     /* We have to check that tng is set because
 173      * tng_util_trajectory_close wants to return a NULL in it, and
 174      * gives a fatal error if it is NULL. */
 175 #ifdef GMX_USE_TNG
 176     if (tng)
 177     {
 178         tng_util_trajectory_close(tng);
 179     }
 180 #else
 181     GMX_UNUSED_VALUE(tng);
 182 #endif
 183 }
 184
 185 #ifdef GMX_USE_TNG
 186 static void addTngMoleculeFromTopology(tng_trajectory_t     tng,
 187                                        const char          *moleculeName,
 188                                        const t_atoms       *atoms,
 189                                        gmx_int64_t          numMolecules,
 190                                        tng_molecule_t      *tngMol)
 191 {
 192     if (tng_molecule_add(tng, moleculeName, tngMol) != TNG_SUCCESS)
 193     {
 194         gmx_file("Cannot add molecule to TNG molecular system.");
 195     }
 196
 197     /* FIXME: The TNG atoms should contain mass and atomB info (for free
 198      * energy calculations), i.e. in when it's available in TNG (2.0). */
 199     for (int atomIt = 0; atomIt < atoms->nr; atomIt++)
 200     {
 201         const t_atom *at = &atoms->atom[atomIt];
 202         /* FIXME: Currently the TNG API can only add atoms belonging to a
 203          * residue and chain. Wait for TNG 2.0*/
 204         if (atoms->nres > 0)
 205         {
 206             const t_resinfo *resInfo        = &atoms->resinfo[at->resind];
 207             char             chainName[2]   = {resInfo->chainid, 0};
 208             tng_chain_t      tngChain       = NULL;
 209             tng_residue_t    tngRes         = NULL;
 210             tng_atom_t       tngAtom        = NULL;
 211
 212             if (tng_molecule_chain_find (tng, *tngMol, chainName,
 213                                          (gmx_int64_t)-1, &tngChain) !=
 214                 TNG_SUCCESS)
 215             {
 216                 tng_molecule_chain_add (tng, *tngMol, chainName,
 217                                         &tngChain);
 218             }
 219
 220             /* FIXME: When TNG supports both residue index and residue
 221              * number the latter should be used. Wait for TNG 2.0*/
 222             if (tng_chain_residue_find(tng, tngChain, *resInfo->name,
 223                                        at->resind + 1, &tngRes)
 224                 != TNG_SUCCESS)
 225             {
 226                 tng_chain_residue_add(tng, tngChain, *resInfo->name, &tngRes);
 227             }
 228             tng_residue_atom_add(tng, tngRes, *(atoms->atomname[atomIt]), *(atoms->atomtype[atomIt]), &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 molIt = 0; molIt < mtop->nmolblock; molIt++)
 248     {
 249         tng_molecule_t       tngMol  = NULL;
 250         const gmx_moltype_t *molType =
 251             &mtop->moltype[mtop->molblock[molIt].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[molIt].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 #ifdef 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         vout        = 0;
 385         fout        = 0;
 386         compression = TNG_TNG_COMPRESSION;
 387     }
 388     else
 389     {
 390         xout        = ir->nstxout;
 391         vout        = ir->nstvout;
 392         fout        = ir->nstfout;
 393         compression = TNG_GZIP_COMPRESSION;
 394     }
 395     if (xout)
 396     {
 397         set_writing_interval(tng, xout, 3, TNG_TRAJ_POSITIONS,
 398                              "POSITIONS", TNG_PARTICLE_BLOCK_DATA,
 399                              compression);
 400         /* The design of TNG makes it awkward to try to write a box
 401          * with multiple periodicities, which might be co-prime. Since
 402          * the use cases for the box with a frame consisting only of
 403          * velocities seem low, for now we associate box writing with
 404          * position writing. */
 405         set_writing_interval(tng, xout, 9, TNG_TRAJ_BOX_SHAPE,
 406                              "BOX SHAPE", TNG_NON_PARTICLE_BLOCK_DATA,
 407                              TNG_GZIP_COMPRESSION);
 408         /* TODO: if/when we write energies to TNG also, reconsider how
 409          * and when box information is written, because GROMACS
 410          * behaviour pre-5.0 was to write the box with every
 411          * trajectory frame and every energy frame, and probably
 412          * people depend on this. */
 413
 414         /* TODO: If we need to write lambda values at steps when
 415          * positions (or other data) are not also being written, then
 416          * code in mdoutf.c will need to match however that is
 417          * organized here. */
 418         set_writing_interval(tng, xout, 1, TNG_GMX_LAMBDA,
 419                              "LAMBDAS", TNG_NON_PARTICLE_BLOCK_DATA,
 420                              TNG_GZIP_COMPRESSION);
 421
 422         /* FIXME: gcd and lowest currently not used. */
 423 //         gcd = greatest_common_divisor_if_positive(gcd, xout);
 424 //         if (lowest < 0 || xout < lowest)
 425 //         {
 426 //             lowest = xout;
 427 //         }
 428     }
 429     if (vout)
 430     {
 431         set_writing_interval(tng, ir->nstvout, 3, TNG_TRAJ_VELOCITIES,
 432                              "VELOCITIES", TNG_PARTICLE_BLOCK_DATA,
 433                              compression);
 434
 435         /* FIXME: gcd and lowest currently not used. */
 436 //         gcd = greatest_common_divisor_if_positive(gcd, vout);
 437 //         if (lowest < 0 || vout < lowest)
 438 //         {
 439 //             lowest = vout;
 440 //         }
 441     }
 442     if (fout)
 443     {
 444         set_writing_interval(tng, ir->nstfout, 3, TNG_TRAJ_FORCES,
 445                              "FORCES", TNG_PARTICLE_BLOCK_DATA,
 446                              TNG_GZIP_COMPRESSION);
 447
 448         /* FIXME: gcd and lowest currently not used. */
 449 //         gcd = greatest_common_divisor_if_positive(gcd, fout);
 450 //         if (lowest < 0 || fout < lowest)
 451 //         {
 452 //             lowest = fout;
 453 //         }
 454     }
 455     /* FIXME: See above. gcd interval for lambdas is disabled. */
 456 //     if (gcd > 0)
 457 //     {
 458 //         /* Lambdas written at an interval of the lowest common denominator
 459 //          * of other output */
 460 //         set_writing_interval(tng, gcd, 1, TNG_GMX_LAMBDA,
 461 //                                  "LAMBDAS", TNG_NON_PARTICLE_BLOCK_DATA,
 462 //                                  TNG_GZIP_COMPRESSION);
 463 //
 464 //         if (gcd < lowest / 10)
 465 //         {
 466 //             gmx_warning("The lowest common denominator of trajectory output is "
 467 //                         "every %d step(s), whereas the shortest output interval "
 468 //                         "is every %d steps.", gcd, lowest);
 469 //         }
 470 //     }
 471 }
 472 #endif
 473
 474 void gmx_tng_prepare_md_writing(tng_trajectory_t  tng,
 475                                 const gmx_mtop_t *mtop,
 476                                 const t_inputrec *ir)
 477 {
 478 #ifdef GMX_USE_TNG
 479     gmx_tng_add_mtop(tng, mtop);
 480     set_writing_intervals(tng, FALSE, ir);
 481     tng_time_per_frame_set(tng, ir->delta_t * PICO);
 482 #else
 483     GMX_UNUSED_VALUE(tng);
 484     GMX_UNUSED_VALUE(mtop);
 485     GMX_UNUSED_VALUE(ir);
 486 #endif
 487 }
 488
 489 #ifdef GMX_USE_TNG
 490 /* Create a TNG molecule representing the selection groups
 491  * to write */
 492 static void add_selection_groups(tng_trajectory_t  tng,
 493                                  const gmx_mtop_t *mtop)
 494 {
 495     const gmx_moltype_t     *molType;
 496     const t_atoms           *atoms;
 497     const t_atom            *at;
 498     const t_resinfo         *resInfo;
 499     const t_ilist           *ilist;
 500     int                      nAtoms      = 0, i = 0, j, molIt, atomIt, nameIndex;
 501     int                      atom_offset = 0;
 502     tng_molecule_t           mol, iterMol;
 503     tng_chain_t              chain;
 504     tng_residue_t            res;
 505     tng_atom_t               atom;
 506     tng_bond_t               tngBond;
 507     gmx_int64_t              nMols;
 508     char                    *groupName;
 509
 510     /* The name of the TNG molecule containing the selection group is the
 511      * same as the name of the selection group. */
 512     nameIndex = *mtop->groups.grps[egcCompressedX].nm_ind;
 513     groupName = *mtop->groups.grpname[nameIndex];
 514
 515     tng_molecule_alloc(tng, &mol);
 516     tng_molecule_name_set(tng, mol, groupName);
 517     tng_molecule_chain_add(tng, mol, "", &chain);
 518     for (molIt = 0; molIt < mtop->nmoltype; molIt++)
 519     {
 520         molType = &mtop->moltype[mtop->molblock[molIt].type];
 521
 522         atoms = &molType->atoms;
 523
 524         for (j = 0; j < mtop->molblock[molIt].nmol; j++)
 525         {
 526             bool bAtomsAdded = FALSE;
 527             for (atomIt = 0; atomIt < atoms->nr; atomIt++, i++)
 528             {
 529                 char *res_name;
 530                 int   res_id;
 531
 532                 if (ggrpnr(&mtop->groups, egcCompressedX, i) != 0)
 533                 {
 534                     continue;
 535                 }
 536                 at = &atoms->atom[atomIt];
 537                 if (atoms->nres > 0)
 538                 {
 539                     resInfo = &atoms->resinfo[at->resind];
 540                     /* FIXME: When TNG supports both residue index and residue
 541                      * number the latter should be used. */
 542                     res_name = *resInfo->name;
 543                     res_id   = at->resind + 1;
 544                 }
 545                 else
 546                 {
 547                     res_name = (char *)"";
 548                     res_id   = 0;
 549                 }
 550                 if (tng_chain_residue_find(tng, chain, res_name, res_id, &res)
 551                     != TNG_SUCCESS)
 552                 {
 553                     /* Since there is ONE chain for selection groups do not keep the
 554                      * original residue IDs - otherwise there might be conflicts. */
 555                     tng_chain_residue_add(tng, chain, res_name, &res);
 556                 }
 557                 tng_residue_atom_w_id_add(tng, res, *(atoms->atomname[atomIt]),
 558                                           *(atoms->atomtype[atomIt]),
 559                                           atom_offset + atomIt, &atom);
 560                 nAtoms++;
 561                 bAtomsAdded = TRUE;
 562             }
 563             /* Add bonds. */
 564             if (bAtomsAdded)
 565             {
 566                 for (int k = 0; k < F_NRE; k++)
 567                 {
 568                     if (IS_CHEMBOND(k))
 569                     {
 570                         ilist = &molType->ilist[k];
 571                         if (ilist)
 572                         {
 573                             int l = 1;
 574                             while (l < ilist->nr)
 575                             {
 576                                 int atom1, atom2;
 577                                 atom1 = ilist->iatoms[l] + atom_offset;
 578                                 atom2 = ilist->iatoms[l+1] + atom_offset;
 579                                 if (ggrpnr(&mtop->groups, egcCompressedX, atom1) == 0 &&
 580                                     ggrpnr(&mtop->groups, egcCompressedX, atom2) == 0)
 581                                 {
 582                                     tng_molecule_bond_add(tng, mol, ilist->iatoms[l],
 583                                                           ilist->iatoms[l+1], &tngBond);
 584                                 }
 585                                 l += 3;
 586                             }
 587                         }
 588                     }
 589                 }
 590                 /* Settle is described using three atoms */
 591                 ilist = &molType->ilist[F_SETTLE];
 592                 if (ilist)
 593                 {
 594                     int l = 1;
 595                     while (l < ilist->nr)
 596                     {
 597                         int atom1, atom2, atom3;
 598                         atom1 = ilist->iatoms[l] + atom_offset;
 599                         atom2 = ilist->iatoms[l+1] + atom_offset;
 600                         atom3 = ilist->iatoms[l+2] + atom_offset;
 601                         if (ggrpnr(&mtop->groups, egcCompressedX, atom1) == 0)
 602                         {
 603                             if (ggrpnr(&mtop->groups, egcCompressedX, atom2) == 0)
 604                             {
 605                                 tng_molecule_bond_add(tng, mol, atom1,
 606                                                       atom2, &tngBond);
 607                             }
 608                             if (ggrpnr(&mtop->groups, egcCompressedX, atom3) == 0)
 609                             {
 610                                 tng_molecule_bond_add(tng, mol, atom1,
 611                                                       atom3, &tngBond);
 612                             }
 613                         }
 614                         l += 4;
 615                     }
 616                 }
 617             }
 618             atom_offset += atoms->nr;
 619         }
 620     }
 621     if (nAtoms != i)
 622     {
 623         tng_molecule_existing_add(tng, &mol);
 624         tng_molecule_cnt_set(tng, mol, 1);
 625         tng_num_molecule_types_get(tng, &nMols);
 626         for (gmx_int64_t k = 0; k < nMols; k++)
 627         {
 628             tng_molecule_of_index_get(tng, k, &iterMol);
 629             if (iterMol == mol)
 630             {
 631                 continue;
 632             }
 633             tng_molecule_cnt_set(tng, iterMol, 0);
 634         }
 635     }
 636     else
 637     {
 638         tng_molecule_free(tng, &mol);
 639     }
 640 }
 641 #endif
 642
 643 void gmx_tng_set_compression_precision(tng_trajectory_t tng,
 644                                        real             prec)
 645 {
 646 #ifdef GMX_USE_TNG
 647     tng_compression_precision_set(tng, 1.0/prec);
 648 #else
 649     GMX_UNUSED_VALUE(tng);
 650     GMX_UNUSED_VALUE(prec);
 651 #endif
 652 }
 653
 654 void gmx_tng_prepare_low_prec_writing(tng_trajectory_t  tng,
 655                                       const gmx_mtop_t *mtop,
 656                                       const t_inputrec *ir)
 657 {
 658 #ifdef GMX_USE_TNG
 659     gmx_tng_add_mtop(tng, mtop);
 660     add_selection_groups(tng, mtop);
 661     set_writing_intervals(tng, TRUE, ir);
 662     tng_time_per_frame_set(tng, ir->delta_t * PICO);
 663     gmx_tng_set_compression_precision(tng, ir->x_compression_precision);
 664 #else
 665     GMX_UNUSED_VALUE(tng);
 666     GMX_UNUSED_VALUE(mtop);
 667     GMX_UNUSED_VALUE(ir);
 668 #endif
 669 }
 670
 671 void gmx_fwrite_tng(tng_trajectory_t tng,
 672                     const gmx_bool   bUseLossyCompression,
 673                     int              step,
 674                     real             elapsedPicoSeconds,
 675                     real             lambda,
 676                     const rvec      *box,
 677                     int              nAtoms,
 678                     const rvec      *x,
 679                     const rvec      *v,
 680                     const rvec      *f)
 681 {
 682 #ifdef GMX_USE_TNG
 683     typedef tng_function_status (*write_data_func_pointer)(tng_trajectory_t,
 684                                                            const gmx_int64_t,
 685                                                            const double,
 686                                                            const real*,
 687                                                            const gmx_int64_t,
 688                                                            const gmx_int64_t,
 689                                                            const char*,
 690                                                            const char,
 691                                                            const char);
 692 #ifdef GMX_DOUBLE
 693     static write_data_func_pointer           write_data           = tng_util_generic_with_time_double_write;
 694 #else
 695     static write_data_func_pointer           write_data           = tng_util_generic_with_time_write;
 696 #endif
 697     double                                   elapsedSeconds = elapsedPicoSeconds * PICO;
 698     gmx_int64_t                              nParticles;
 699     char                                     compression;
 700
 701
 702     if (!tng)
 703     {
 704         /* This function might get called when the type of the
 705            compressed trajectory is actually XTC. So we exit and move
 706            on. */
 707         return;
 708     }
 709
 710     tng_num_particles_get(tng, &nParticles);
 711     if (nAtoms != (int)nParticles)
 712     {
 713         tng_implicit_num_particles_set(tng, nAtoms);
 714     }
 715
 716     if (bUseLossyCompression)
 717     {
 718         compression = TNG_TNG_COMPRESSION;
 719     }
 720     else
 721     {
 722         compression = TNG_GZIP_COMPRESSION;
 723     }
 724
 725     /* The writing is done using write_data, which writes float or double
 726      * depending on the GROMACS compilation. */
 727     if (x)
 728     {
 729         GMX_ASSERT(box, "Need a non-NULL box if positions are written");
 730
 731         if (write_data(tng, step, elapsedSeconds,
 732                        reinterpret_cast<const real *>(x),
 733                        3, TNG_TRAJ_POSITIONS, "POSITIONS",
 734                        TNG_PARTICLE_BLOCK_DATA,
 735                        compression) != TNG_SUCCESS)
 736         {
 737             gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 738         }
 739         /* TNG-MF1 compression only compresses positions and velocities. Use lossless
 740          * compression for box shape regardless of output mode */
 741         if (write_data(tng, step, elapsedSeconds,
 742                        reinterpret_cast<const real *>(box),
 743                        9, TNG_TRAJ_BOX_SHAPE, "BOX SHAPE",
 744                        TNG_NON_PARTICLE_BLOCK_DATA,
 745                        TNG_GZIP_COMPRESSION) != TNG_SUCCESS)
 746         {
 747             gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 748         }
 749     }
 750
 751     if (v)
 752     {
 753         if (write_data(tng, step, elapsedSeconds,
 754                        reinterpret_cast<const real *>(v),
 755                        3, TNG_TRAJ_VELOCITIES, "VELOCITIES",
 756                        TNG_PARTICLE_BLOCK_DATA,
 757                        compression) != TNG_SUCCESS)
 758         {
 759             gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 760         }
 761     }
 762
 763     if (f)
 764     {
 765         /* TNG-MF1 compression only compresses positions and velocities. Use lossless
 766          * compression for forces regardless of output mode */
 767         if (write_data(tng, step, elapsedSeconds,
 768                        reinterpret_cast<const real *>(f),
 769                        3, TNG_TRAJ_FORCES, "FORCES",
 770                        TNG_PARTICLE_BLOCK_DATA,
 771                        TNG_GZIP_COMPRESSION) != TNG_SUCCESS)
 772         {
 773             gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 774         }
 775     }
 776
 777     /* TNG-MF1 compression only compresses positions and velocities. Use lossless
 778      * compression for lambdas regardless of output mode */
 779     if (write_data(tng, step, elapsedSeconds,
 780                    reinterpret_cast<const real *>(&lambda),
 781                    1, TNG_GMX_LAMBDA, "LAMBDAS",
 782                    TNG_NON_PARTICLE_BLOCK_DATA,
 783                    TNG_GZIP_COMPRESSION) != TNG_SUCCESS)
 784     {
 785         gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
 786     }
 787 #else
 788     GMX_UNUSED_VALUE(tng);
 789     GMX_UNUSED_VALUE(bUseLossyCompression);
 790     GMX_UNUSED_VALUE(step);
 791     GMX_UNUSED_VALUE(elapsedPicoSeconds);
 792     GMX_UNUSED_VALUE(lambda);
 793     GMX_UNUSED_VALUE(box);
 794     GMX_UNUSED_VALUE(nAtoms);
 795     GMX_UNUSED_VALUE(x);
 796     GMX_UNUSED_VALUE(v);
 797     GMX_UNUSED_VALUE(f);
 798 #endif
 799 }
 800
 801 void fflush_tng(tng_trajectory_t tng)
 802 {
 803 #ifdef GMX_USE_TNG
 804     if (!tng)
 805     {
 806         return;
 807     }
 808     tng_frame_set_premature_write(tng, TNG_USE_HASH);
 809 #else
 810     GMX_UNUSED_VALUE(tng);
 811 #endif
 812 }
 813
 814 float gmx_tng_get_time_of_final_frame(tng_trajectory_t tng)
 815 {
 816 #ifdef GMX_USE_TNG
 817     gmx_int64_t nFrames;
 818     double      time;
 819     float       fTime;
 820
 821     tng_num_frames_get(tng, &nFrames);
 822     tng_util_time_of_frame_get(tng, nFrames - 1, &time);
 823
 824     fTime = time / PICO;
 825     return fTime;
 826 #else
 827     GMX_UNUSED_VALUE(tng);
 828     return -1.0;
 829 #endif
 830 }