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