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