src/gromacs/fileio/tngio.cpp

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