src/gromacs/fileio/tngio.cpp

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