1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 2013,2014,2015,2018,2019,2020, 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 #include <cmath>
42
43 #include <algorithm>
44 #include <iterator>
45 #include <memory>
46 #include <numeric>
47 #include <vector>
48
49 #if GMX_USE_TNG
50 # include "tng/tng_io.h"
51 #endif
52
53 #include "gromacs/math/units.h"
54 #include "gromacs/math/utilities.h"
55 #include "gromacs/mdtypes/inputrec.h"
56 #include "gromacs/topology/ifunc.h"
57 #include "gromacs/topology/topology.h"
58 #include "gromacs/trajectory/trajectoryframe.h"
59 #include "gromacs/utility/basedefinitions.h"
60 #include "gromacs/utility/baseversion.h"
61 #include "gromacs/utility/fatalerror.h"
62 #include "gromacs/utility/futil.h"
63 #include "gromacs/utility/gmxassert.h"
64 #include "gromacs/utility/programcontext.h"
65 #include "gromacs/utility/smalloc.h"
66 #include "gromacs/utility/sysinfo.h"
67 #include "gromacs/utility/unique_cptr.h"
68
69 #if !GMX_USE_TNG
70 using tng_trajectory_t = void*;
71 #endif
72
73 /*! \brief Gromacs Wrapper around tng datatype
74 *
75 * This could in principle hold any GROMACS-specific requirements not yet
76 * implemented in or not relevant to the TNG library itself. However, for now
77 * we only use it to handle some shortcomings we have discovered, where the TNG
78 * API itself is a bit fragile and can end up overwriting data if called several
79 * times with the same frame number.
80 * The logic to determine the time per step was also a bit fragile. This is not
81 * critical, but since we anyway need a wrapper for ensuring unique frame
82 * numbers, we can also use it to store the time of the first step and use that
83 * to derive a slightly better/safer estimate of the time per step.
84 *
85 * At some future point where we have a second-generation TNG API we should
86 * consider removing this again.
87 */
88 struct gmx_tng_trajectory
89 {
90 tng_trajectory_t tng; //!< Actual TNG handle (pointer)
91 bool lastStepDataIsValid; //!< True if lastStep has been set
92 std::int64_t lastStep; //!< Index/step used for last frame
93 bool lastTimeDataIsValid; //!< True if lastTime has been set
94 double lastTime; //!< Time of last frame (TNG unit is seconds)
95 bool timePerFrameIsSet; //!< True if we have set the time per frame
96 int boxOutputInterval; //!< Number of steps between the output of box size
97 int lambdaOutputInterval; //!< Number of steps between the output of lambdas
98 };
99
100 #if GMX_USE_TNG
modeToVerb(char mode)101 static const char* modeToVerb(char mode)
102 {
103 const char* p;
104 switch (mode)
105 {
106 case 'r': p = "reading"; break;
107 case 'w': p = "writing"; break;
108 case 'a': p = "appending"; break;
109 default: gmx_fatal(FARGS, "Invalid file opening mode %c", mode);
110 }
111 return p;
112 }
113 #endif
114
gmx_tng_open(const char * filename,char mode,gmx_tng_trajectory_t * gmx_tng)115 void gmx_tng_open(const char* filename, char mode, gmx_tng_trajectory_t* gmx_tng)
116 {
117 #if GMX_USE_TNG
118 /* First check whether we have to make a backup,
119 * only for writing, not for read or append.
120 */
121 if (mode == 'w')
122 {
123 make_backup(filename);
124 }
125
126 *gmx_tng = new gmx_tng_trajectory;
127 (*gmx_tng)->lastStepDataIsValid = false;
128 (*gmx_tng)->lastTimeDataIsValid = false;
129 (*gmx_tng)->timePerFrameIsSet = false;
130 tng_trajectory_t* tng = &(*gmx_tng)->tng;
131
132 /* tng must not be pointing at already allocated memory.
133 * Memory will be allocated by tng_util_trajectory_open() and must
134 * later on be freed by tng_util_trajectory_close(). */
135 if (TNG_SUCCESS != tng_util_trajectory_open(filename, mode, tng))
136 {
137 /* TNG does return more than one degree of error, but there is
138 no use case for GROMACS handling the non-fatal errors
139 gracefully. */
140 gmx_fatal(FARGS, "File I/O error while opening %s for %s", filename, modeToVerb(mode));
141 }
142
143 if (mode == 'w' || mode == 'a')
144 {
145 char hostname[256];
146 gmx_gethostname(hostname, 256);
147 if (mode == 'w')
148 {
149 tng_first_computer_name_set(*tng, hostname);
150 }
151 else
152 {
153 tng_last_computer_name_set(*tng, hostname);
154 }
155
156 char programInfo[256];
157 const char* precisionString = "";
158 # if GMX_DOUBLE
159 precisionString = " (double precision)";
160 # endif
161 sprintf(programInfo, "%.100s %.128s%.24s", gmx::getProgramContext().displayName(),
162 gmx_version(), precisionString);
163 if (mode == 'w')
164 {
165 tng_first_program_name_set(*tng, programInfo);
166 }
167 else
168 {
169 tng_last_program_name_set(*tng, programInfo);
170 }
171
172 char username[256];
173 if (!gmx_getusername(username, 256))
174 {
175 if (mode == 'w')
176 {
177 tng_first_user_name_set(*tng, username);
178 }
179 else
180 {
181 tng_last_user_name_set(*tng, username);
182 tng_file_headers_write(*tng, TNG_USE_HASH);
183 }
184 }
185 }
186 #else
187 gmx_file("GROMACS was compiled without TNG support, cannot handle this file type");
188 GMX_UNUSED_VALUE(filename);
189 GMX_UNUSED_VALUE(mode);
190 GMX_UNUSED_VALUE(gmx_tng);
191 #endif
192 }
193
gmx_tng_close(gmx_tng_trajectory_t * gmx_tng)194 void gmx_tng_close(gmx_tng_trajectory_t* gmx_tng)
195 {
196 /* We have to check that tng is set because
197 * tng_util_trajectory_close wants to return a NULL in it, and
198 * gives a fatal error if it is NULL. */
199 #if GMX_USE_TNG
200 if (gmx_tng == nullptr || *gmx_tng == nullptr)
201 {
202 return;
203 }
204 tng_trajectory_t* tng = &(*gmx_tng)->tng;
205
206 if (tng)
207 {
208 tng_util_trajectory_close(tng);
209 }
210 delete *gmx_tng;
211 *gmx_tng = nullptr;
212
213 #else
214 GMX_UNUSED_VALUE(gmx_tng);
215 #endif
216 }
217
218 #if GMX_USE_TNG
addTngMoleculeFromTopology(gmx_tng_trajectory_t gmx_tng,const char * moleculeName,const t_atoms * atoms,int64_t numMolecules,tng_molecule_t * tngMol)219 static void addTngMoleculeFromTopology(gmx_tng_trajectory_t gmx_tng,
220 const char* moleculeName,
221 const t_atoms* atoms,
222 int64_t numMolecules,
223 tng_molecule_t* tngMol)
224 {
225 tng_trajectory_t tng = gmx_tng->tng;
226 tng_chain_t tngChain = nullptr;
227 tng_residue_t tngRes = nullptr;
228
229 if (tng_molecule_add(tng, moleculeName, tngMol) != TNG_SUCCESS)
230 {
231 gmx_file("Cannot add molecule to TNG molecular system.");
232 }
233
234 for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++)
235 {
236 const t_atom* at = &atoms->atom[atomIndex];
237 /* FIXME: Currently the TNG API can only add atoms belonging to a
238 * residue and chain. Wait for TNG 2.0*/
239 if (atoms->nres > 0)
240 {
241 const t_resinfo* resInfo = &atoms->resinfo[at->resind];
242 char chainName[2] = { resInfo->chainid, 0 };
243 tng_atom_t tngAtom = nullptr;
244 t_atom* prevAtom;
245
246 if (atomIndex > 0)
247 {
248 prevAtom = &atoms->atom[atomIndex - 1];
249 }
250 else
251 {
252 prevAtom = nullptr;
253 }
254
255 /* If this is the first atom or if the residue changed add the
256 * residue to the TNG molecular system. */
257 if (!prevAtom || resInfo != &atoms->resinfo[prevAtom->resind])
258 {
259 /* If this is the first atom or if the chain changed add
260 * the chain to the TNG molecular system. */
261 if (!prevAtom || resInfo->chainid != atoms->resinfo[prevAtom->resind].chainid)
262 {
263 tng_molecule_chain_add(tng, *tngMol, chainName, &tngChain);
264 }
265 /* FIXME: When TNG supports both residue index and residue
266 * number the latter should be used. Wait for TNG 2.0*/
267 tng_chain_residue_add(tng, tngChain, *resInfo->name, &tngRes);
268 }
269 tng_residue_atom_add(tng, tngRes, *(atoms->atomname[atomIndex]),
270 *(atoms->atomtype[atomIndex]), &tngAtom);
271 }
272 }
273 tng_molecule_cnt_set(tng, *tngMol, numMolecules);
274 }
275
gmx_tng_add_mtop(gmx_tng_trajectory_t gmx_tng,const gmx_mtop_t * mtop)276 void gmx_tng_add_mtop(gmx_tng_trajectory_t gmx_tng, const gmx_mtop_t* mtop)
277 {
278 int i;
279 int j;
280 std::vector<real> atomCharges;
281 std::vector<real> atomMasses;
282 tng_bond_t tngBond;
283 char datatype;
284
285 tng_trajectory_t tng = gmx_tng->tng;
286
287 if (!mtop)
288 {
289 /* No topology information available to add. */
290 return;
291 }
292
293 # if GMX_DOUBLE
294 datatype = TNG_DOUBLE_DATA;
295 # else
296 datatype = TNG_FLOAT_DATA;
297 # endif
298
299 atomCharges.reserve(mtop->natoms);
300 atomMasses.reserve(mtop->natoms);
301
302 for (const gmx_molblock_t& molBlock : mtop->molblock)
303 {
304 tng_molecule_t tngMol = nullptr;
305 const gmx_moltype_t* molType = &mtop->moltype[molBlock.type];
306
307 /* Add a molecule to the TNG trajectory with the same name as the
308 * current molecule. */
309 addTngMoleculeFromTopology(gmx_tng, *(molType->name), &molType->atoms, molBlock.nmol, &tngMol);
310
311 /* Bonds have to be deduced from interactions (constraints etc). Different
312 * interactions have different sets of parameters. */
313 /* Constraints are specified using two atoms */
314 for (i = 0; i < F_NRE; i++)
315 {
316 if (IS_CHEMBOND(i))
317 {
318 const InteractionList& ilist = molType->ilist[i];
319 j = 1;
320 while (j < ilist.size())
321 {
322 tng_molecule_bond_add(tng, tngMol, ilist.iatoms[j], ilist.iatoms[j + 1], &tngBond);
323 j += 3;
324 }
325 }
326 }
327 /* Settle is described using three atoms */
328 const InteractionList& ilist = molType->ilist[F_SETTLE];
329 j = 1;
330 while (j < ilist.size())
331 {
332 tng_molecule_bond_add(tng, tngMol, ilist.iatoms[j], ilist.iatoms[j + 1], &tngBond);
333 tng_molecule_bond_add(tng, tngMol, ilist.iatoms[j], ilist.iatoms[j + 2], &tngBond);
334 j += 4;
335 }
336 /* First copy atom charges and masses, first atom by atom and then copy the memory for the molecule instances.
337 * FIXME: Atom B state data should also be written to TNG (v 2.0?) */
338 for (int atomCounter = 0; atomCounter < molType->atoms.nr; atomCounter++)
339 {
340 atomCharges.push_back(molType->atoms.atom[atomCounter].q);
341 atomMasses.push_back(molType->atoms.atom[atomCounter].m);
342 }
343 for (int molCounter = 1; molCounter < molBlock.nmol; molCounter++)
344 {
345 std::copy_n(atomCharges.end() - molType->atoms.nr, molType->atoms.nr,
346 std::back_inserter(atomCharges));
347 std::copy_n(atomMasses.end() - molType->atoms.nr, molType->atoms.nr,
348 std::back_inserter(atomMasses));
349 }
350 }
351 /* Write the TNG data blocks. */
352 tng_particle_data_block_add(tng, TNG_TRAJ_PARTIAL_CHARGES, "PARTIAL CHARGES", datatype,
353 TNG_NON_TRAJECTORY_BLOCK, 1, 1, 1, 0, mtop->natoms,
354 TNG_GZIP_COMPRESSION, atomCharges.data());
355 tng_particle_data_block_add(tng, TNG_TRAJ_MASSES, "ATOM MASSES", datatype, TNG_NON_TRAJECTORY_BLOCK,
356 1, 1, 1, 0, mtop->natoms, TNG_GZIP_COMPRESSION, atomMasses.data());
357 }
358
359 /*! \libinternal \brief Compute greatest common divisor of n1 and n2
360 * if they are positive.
361 *
362 * If only one of n1 and n2 is positive, then return it.
363 * If neither n1 or n2 is positive, then return -1. */
greatest_common_divisor_if_positive(int n1,int n2)364 static int greatest_common_divisor_if_positive(int n1, int n2)
365 {
366 if (0 >= n1)
367 {
368 return (0 >= n2) ? -1 : n2;
369 }
370 if (0 >= n2)
371 {
372 return n1;
373 }
374
375 /* We have a non-trivial greatest common divisor to compute. */
376 return std::gcd(n1, n2);
377 }
378
379 /* By default try to write 100 frames (of actual output) in each frame set.
380 * This number is the number of outputs of the most frequently written data
381 * type per frame set.
382 * TODO for 5.1: Verify that 100 frames per frame set is efficient for most
383 * setups regarding compression efficiency and compression time. Make this
384 * a hidden command-line option? */
385 const int defaultFramesPerFrameSet = 100;
386
387 /*! \libinternal \brief Set the number of frames per frame
388 * set according to output intervals.
389 * The default is that 100 frames are written of the data
390 * that is written most often. */
tng_set_frames_per_frame_set(gmx_tng_trajectory_t gmx_tng,const gmx_bool bUseLossyCompression,const t_inputrec * ir)391 static void tng_set_frames_per_frame_set(gmx_tng_trajectory_t gmx_tng,
392 const gmx_bool bUseLossyCompression,
393 const t_inputrec* ir)
394 {
395 int gcd = -1;
396 tng_trajectory_t tng = gmx_tng->tng;
397
398 /* Set the number of frames per frame set to contain at least
399 * defaultFramesPerFrameSet of the lowest common denominator of
400 * the writing interval of positions and velocities. */
401 /* FIXME after 5.0: consider nstenergy also? */
402 if (bUseLossyCompression)
403 {
404 gcd = ir->nstxout_compressed;
405 }
406 else
407 {
408 gcd = greatest_common_divisor_if_positive(ir->nstxout, ir->nstvout);
409 gcd = greatest_common_divisor_if_positive(gcd, ir->nstfout);
410 }
411 if (0 >= gcd)
412 {
413 return;
414 }
415
416 tng_num_frames_per_frame_set_set(tng, gcd * defaultFramesPerFrameSet);
417 }
418
419 /*! \libinternal \brief Set the data-writing intervals, and number of
420 * frames per frame set */
set_writing_intervals(gmx_tng_trajectory_t gmx_tng,const gmx_bool bUseLossyCompression,const t_inputrec * ir)421 static void set_writing_intervals(gmx_tng_trajectory_t gmx_tng,
422 const gmx_bool bUseLossyCompression,
423 const t_inputrec* ir)
424 {
425 tng_trajectory_t tng = gmx_tng->tng;
426
427 /* Define pointers to specific writing functions depending on if we
428 * write float or double data */
429 typedef tng_function_status (*set_writing_interval_func_pointer)(
430 tng_trajectory_t, const int64_t, const int64_t, const int64_t, const char*, const char,
431 const char);
432 # if GMX_DOUBLE
433 set_writing_interval_func_pointer set_writing_interval = tng_util_generic_write_interval_double_set;
434 # else
435 set_writing_interval_func_pointer set_writing_interval = tng_util_generic_write_interval_set;
436 # endif
437 int xout, vout, fout;
438 int gcd = -1, lowest = -1;
439 char compression;
440
441 tng_set_frames_per_frame_set(gmx_tng, bUseLossyCompression, ir);
442
443 if (bUseLossyCompression)
444 {
445 xout = ir->nstxout_compressed;
446
447 /* If there is no uncompressed coordinate output write forces
448 and velocities to the compressed tng file. */
449 if (ir->nstxout)
450 {
451 vout = 0;
452 fout = 0;
453 }
454 else
455 {
456 vout = ir->nstvout;
457 fout = ir->nstfout;
458 }
459 compression = TNG_TNG_COMPRESSION;
460 }
461 else
462 {
463 xout = ir->nstxout;
464 vout = ir->nstvout;
465 fout = ir->nstfout;
466 compression = TNG_GZIP_COMPRESSION;
467 }
468 if (xout)
469 {
470 set_writing_interval(tng, xout, 3, TNG_TRAJ_POSITIONS, "POSITIONS", TNG_PARTICLE_BLOCK_DATA,
471 compression);
472 /* TODO: if/when we write energies to TNG also, reconsider how
473 * and when box information is written, because GROMACS
474 * behaviour pre-5.0 was to write the box with every
475 * trajectory frame and every energy frame, and probably
476 * people depend on this. */
477
478 gcd = greatest_common_divisor_if_positive(gcd, xout);
479 if (lowest < 0 || xout < lowest)
480 {
481 lowest = xout;
482 }
483 }
484 if (vout)
485 {
486 set_writing_interval(tng, vout, 3, TNG_TRAJ_VELOCITIES, "VELOCITIES",
487 TNG_PARTICLE_BLOCK_DATA, compression);
488
489 gcd = greatest_common_divisor_if_positive(gcd, vout);
490 if (lowest < 0 || vout < lowest)
491 {
492 lowest = vout;
493 }
494 }
495 if (fout)
496 {
497 set_writing_interval(tng, fout, 3, TNG_TRAJ_FORCES, "FORCES", TNG_PARTICLE_BLOCK_DATA,
498 TNG_GZIP_COMPRESSION);
499
500 gcd = greatest_common_divisor_if_positive(gcd, fout);
501 if (lowest < 0 || fout < lowest)
502 {
503 lowest = fout;
504 }
505 }
506 if (gcd > 0)
507 {
508 /* Lambdas and box shape written at an interval of the lowest common
509 denominator of other output */
510 set_writing_interval(tng, gcd, 1, TNG_GMX_LAMBDA, "LAMBDAS", TNG_NON_PARTICLE_BLOCK_DATA,
511 TNG_GZIP_COMPRESSION);
512
513 set_writing_interval(tng, gcd, 9, TNG_TRAJ_BOX_SHAPE, "BOX SHAPE",
514 TNG_NON_PARTICLE_BLOCK_DATA, TNG_GZIP_COMPRESSION);
515 gmx_tng->lambdaOutputInterval = gcd;
516 gmx_tng->boxOutputInterval = gcd;
517 if (gcd < lowest / 10)
518 {
519 gmx_warning(
520 "The lowest common denominator of trajectory output is "
521 "every %d step(s), whereas the shortest output interval "
522 "is every %d steps.",
523 gcd, lowest);
524 }
525 }
526 }
527 #endif
528
gmx_tng_prepare_md_writing(gmx_tng_trajectory_t gmx_tng,const gmx_mtop_t * mtop,const t_inputrec * ir)529 void gmx_tng_prepare_md_writing(gmx_tng_trajectory_t gmx_tng, const gmx_mtop_t* mtop, const t_inputrec* ir)
530 {
531 #if GMX_USE_TNG
532 gmx_tng_add_mtop(gmx_tng, mtop);
533 set_writing_intervals(gmx_tng, FALSE, ir);
534 tng_time_per_frame_set(gmx_tng->tng, ir->delta_t * PICO);
535 gmx_tng->timePerFrameIsSet = true;
536 #else
537 GMX_UNUSED_VALUE(gmx_tng);
538 GMX_UNUSED_VALUE(mtop);
539 GMX_UNUSED_VALUE(ir);
540 #endif
541 }
542
543 #if GMX_USE_TNG
544 /* Check if all atoms in the molecule system are selected
545 * by a selection group of type specified by gtype. */
all_atoms_selected(const gmx_mtop_t * mtop,const SimulationAtomGroupType gtype)546 static gmx_bool all_atoms_selected(const gmx_mtop_t* mtop, const SimulationAtomGroupType gtype)
547 {
548 /* Iterate through all atoms in the molecule system and
549 * check if they belong to a selection group of the
550 * requested type. */
551 int i = 0;
552 for (const gmx_molblock_t& molBlock : mtop->molblock)
553 {
554 const gmx_moltype_t& molType = mtop->moltype[molBlock.type];
555 const t_atoms& atoms = molType.atoms;
556
557 for (int j = 0; j < molBlock.nmol; j++)
558 {
559 for (int atomIndex = 0; atomIndex < atoms.nr; atomIndex++, i++)
560 {
561 if (getGroupType(mtop->groups, gtype, i) != 0)
562 {
563 return FALSE;
564 }
565 }
566 }
567 }
568 return TRUE;
569 }
570
571 /* Create TNG molecules which will represent each of the selection
572 * groups for writing. But do that only if there is actually a
573 * specified selection group and it is not the whole system.
574 * TODO: Currently the only selection that is taken into account
575 * is egcCompressedX, but other selections should be added when
576 * e.g. writing energies is implemented.
577 */
add_selection_groups(gmx_tng_trajectory_t gmx_tng,const gmx_mtop_t * mtop)578 static void add_selection_groups(gmx_tng_trajectory_t gmx_tng, const gmx_mtop_t* mtop)
579 {
580 const t_atoms* atoms;
581 const t_atom* at;
582 const t_resinfo* resInfo;
583 int nameIndex;
584 int atom_offset = 0;
585 tng_molecule_t mol, iterMol;
586 tng_chain_t chain;
587 tng_residue_t res;
588 tng_atom_t atom;
589 tng_bond_t tngBond;
590 int64_t nMols;
591 char* groupName;
592 tng_trajectory_t tng = gmx_tng->tng;
593
594 /* TODO: When the TNG molecules block is more flexible TNG selection
595 * groups should not need all atoms specified. It should be possible
596 * just to specify what molecules are selected (and/or which atoms in
597 * the molecule) and how many (if applicable). */
598
599 /* If no atoms are selected we do not need to create a
600 * TNG selection group molecule. */
601 if (mtop->groups.numberOfGroupNumbers(SimulationAtomGroupType::CompressedPositionOutput) == 0)
602 {
603 return;
604 }
605
606 /* If all atoms are selected we do not have to create a selection
607 * group molecule in the TNG molecule system. */
608 if (all_atoms_selected(mtop, SimulationAtomGroupType::CompressedPositionOutput))
609 {
610 return;
611 }
612
613 /* The name of the TNG molecule containing the selection group is the
614 * same as the name of the selection group. */
615 nameIndex = mtop->groups.groups[SimulationAtomGroupType::CompressedPositionOutput][0];
616 groupName = *mtop->groups.groupNames[nameIndex];
617
618 tng_molecule_alloc(tng, &mol);
619 tng_molecule_name_set(tng, mol, groupName);
620 tng_molecule_chain_add(tng, mol, "", &chain);
621 int i = 0;
622 for (const gmx_molblock_t& molBlock : mtop->molblock)
623 {
624 const gmx_moltype_t& molType = mtop->moltype[molBlock.type];
625
626 atoms = &molType.atoms;
627
628 for (int j = 0; j < molBlock.nmol; j++)
629 {
630 bool bAtomsAdded = FALSE;
631 for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++, i++)
632 {
633 char* res_name;
634 int res_id;
635
636 if (getGroupType(mtop->groups, SimulationAtomGroupType::CompressedPositionOutput, i) != 0)
637 {
638 continue;
639 }
640 at = &atoms->atom[atomIndex];
641 if (atoms->nres > 0)
642 {
643 resInfo = &atoms->resinfo[at->resind];
644 /* FIXME: When TNG supports both residue index and residue
645 * number the latter should be used. */
646 res_name = *resInfo->name;
647 res_id = at->resind + 1;
648 }
649 else
650 {
651 res_name = const_cast<char*>("");
652 res_id = 0;
653 }
654 if (tng_chain_residue_find(tng, chain, res_name, res_id, &res) != TNG_SUCCESS)
655 {
656 /* Since there is ONE chain for selection groups do not keep the
657 * original residue IDs - otherwise there might be conflicts. */
658 tng_chain_residue_add(tng, chain, res_name, &res);
659 }
660 tng_residue_atom_w_id_add(tng, res, *(atoms->atomname[atomIndex]),
661 *(atoms->atomtype[atomIndex]), atom_offset + atomIndex, &atom);
662 bAtomsAdded = TRUE;
663 }
664 /* Add bonds. */
665 if (bAtomsAdded)
666 {
667 for (int k = 0; k < F_NRE; k++)
668 {
669 if (IS_CHEMBOND(k))
670 {
671 const InteractionList& ilist = molType.ilist[k];
672 for (int l = 1; l < ilist.size(); l += 3)
673 {
674 int atom1, atom2;
675 atom1 = ilist.iatoms[l] + atom_offset;
676 atom2 = ilist.iatoms[l + 1] + atom_offset;
677 if (getGroupType(mtop->groups,
678 SimulationAtomGroupType::CompressedPositionOutput, atom1)
679 == 0
680 && getGroupType(mtop->groups,
681 SimulationAtomGroupType::CompressedPositionOutput, atom2)
682 == 0)
683 {
684 tng_molecule_bond_add(tng, mol, ilist.iatoms[l],
685 ilist.iatoms[l + 1], &tngBond);
686 }
687 }
688 }
689 }
690 /* Settle is described using three atoms */
691 const InteractionList& ilist = molType.ilist[F_SETTLE];
692 for (int l = 1; l < ilist.size(); l += 4)
693 {
694 int atom1, atom2, atom3;
695 atom1 = ilist.iatoms[l] + atom_offset;
696 atom2 = ilist.iatoms[l + 1] + atom_offset;
697 atom3 = ilist.iatoms[l + 2] + atom_offset;
698 if (getGroupType(mtop->groups, SimulationAtomGroupType::CompressedPositionOutput, atom1)
699 == 0)
700 {
701 if (getGroupType(mtop->groups, SimulationAtomGroupType::CompressedPositionOutput, atom2)
702 == 0)
703 {
704 tng_molecule_bond_add(tng, mol, atom1, atom2, &tngBond);
705 }
706 if (getGroupType(mtop->groups, SimulationAtomGroupType::CompressedPositionOutput, atom3)
707 == 0)
708 {
709 tng_molecule_bond_add(tng, mol, atom1, atom3, &tngBond);
710 }
711 }
712 }
713 }
714 atom_offset += atoms->nr;
715 }
716 }
717 tng_molecule_existing_add(tng, &mol);
718 tng_molecule_cnt_set(tng, mol, 1);
719 tng_num_molecule_types_get(tng, &nMols);
720 for (int64_t k = 0; k < nMols; k++)
721 {
722 tng_molecule_of_index_get(tng, k, &iterMol);
723 if (iterMol == mol)
724 {
725 continue;
726 }
727 tng_molecule_cnt_set(tng, iterMol, 0);
728 }
729 }
730 #endif
731
gmx_tng_set_compression_precision(gmx_tng_trajectory_t gmx_tng,real prec)732 void gmx_tng_set_compression_precision(gmx_tng_trajectory_t gmx_tng, real prec)
733 {
734 #if GMX_USE_TNG
735 tng_compression_precision_set(gmx_tng->tng, prec);
736 #else
737 GMX_UNUSED_VALUE(gmx_tng);
738 GMX_UNUSED_VALUE(prec);
739 #endif
740 }
741
gmx_tng_prepare_low_prec_writing(gmx_tng_trajectory_t gmx_tng,const gmx_mtop_t * mtop,const t_inputrec * ir)742 void gmx_tng_prepare_low_prec_writing(gmx_tng_trajectory_t gmx_tng, const gmx_mtop_t* mtop, const t_inputrec* ir)
743 {
744 #if GMX_USE_TNG
745 gmx_tng_add_mtop(gmx_tng, mtop);
746 add_selection_groups(gmx_tng, mtop);
747 set_writing_intervals(gmx_tng, TRUE, ir);
748 tng_time_per_frame_set(gmx_tng->tng, ir->delta_t * PICO);
749 gmx_tng->timePerFrameIsSet = true;
750 gmx_tng_set_compression_precision(gmx_tng, ir->x_compression_precision);
751 #else
752 GMX_UNUSED_VALUE(gmx_tng);
753 GMX_UNUSED_VALUE(mtop);
754 GMX_UNUSED_VALUE(ir);
755 #endif
756 }
757
gmx_fwrite_tng(gmx_tng_trajectory_t gmx_tng,const gmx_bool bUseLossyCompression,int64_t step,real elapsedPicoSeconds,real lambda,const rvec * box,int nAtoms,const rvec * x,const rvec * v,const rvec * f)758 void gmx_fwrite_tng(gmx_tng_trajectory_t gmx_tng,
759 const gmx_bool bUseLossyCompression,
760 int64_t step,
761 real elapsedPicoSeconds,
762 real lambda,
763 const rvec* box,
764 int nAtoms,
765 const rvec* x,
766 const rvec* v,
767 const rvec* f)
768 {
769 #if GMX_USE_TNG
770 typedef tng_function_status (*write_data_func_pointer)(
771 tng_trajectory_t, const int64_t, const double, const real*, const int64_t,
772 const int64_t, const char*, const char, const char);
773 # if GMX_DOUBLE
774 static write_data_func_pointer write_data = tng_util_generic_with_time_double_write;
775 # else
776 static write_data_func_pointer write_data = tng_util_generic_with_time_write;
777 # endif
778 double elapsedSeconds = elapsedPicoSeconds * PICO;
779 int64_t nParticles;
780 char compression;
781
782
783 if (!gmx_tng)
784 {
785 /* This function might get called when the type of the
786 compressed trajectory is actually XTC. So we exit and move
787 on. */
788 return;
789 }
790 tng_trajectory_t tng = gmx_tng->tng;
791
792 // While the GROMACS interface to this routine specifies 'step', TNG itself
793 // only uses 'frame index' internally, although it suggests that it's a good
794 // idea to let this represent the step rather than just frame numbers.
795 //
796 // However, the way the GROMACS interface works, there are cases where
797 // the step is simply not set, so all frames rather have step=0.
798 // If we call the lower-level TNG routines with the same frame index
799 // (which is set from the step) multiple times, things will go very
800 // wrong (overwritten data).
801 //
802 // To avoid this, we require the step value to always be larger than the
803 // previous value, and if this is not true we simply set it to a value
804 // one unit larger than the previous step.
805 //
806 // Note: We do allow the user to specify any crazy value the want for the
807 // first step. If it is "not set", GROMACS will have used the value 0.
808 if (gmx_tng->lastStepDataIsValid && step <= gmx_tng->lastStep)
809 {
810 step = gmx_tng->lastStep + 1;
811 }
812
813 // Now that we have fixed the potentially incorrect step, we can also set
814 // the time per frame if it was not already done, and we have
815 // step/time information from the last frame so it is possible to calculate it.
816 // Note that we do allow the time to be the same for all steps, or even
817 // decreasing. In that case we will simply say the time per step is 0
818 // or negative. A bit strange, but a correct representation of the data :-)
819 if (!gmx_tng->timePerFrameIsSet && gmx_tng->lastTimeDataIsValid && gmx_tng->lastStepDataIsValid)
820 {
821 double deltaTime = elapsedSeconds - gmx_tng->lastTime;
822 std::int64_t deltaStep = step - gmx_tng->lastStep;
823 tng_time_per_frame_set(tng, deltaTime / deltaStep);
824 gmx_tng->timePerFrameIsSet = true;
825 }
826
827 tng_num_particles_get(tng, &nParticles);
828 if (nAtoms != static_cast<int>(nParticles))
829 {
830 tng_implicit_num_particles_set(tng, nAtoms);
831 }
832
833 if (bUseLossyCompression)
834 {
835 compression = TNG_TNG_COMPRESSION;
836 }
837 else
838 {
839 compression = TNG_GZIP_COMPRESSION;
840 }
841
842 /* The writing is done using write_data, which writes float or double
843 * depending on the GROMACS compilation. */
844 if (x)
845 {
846 GMX_ASSERT(box, "Need a non-NULL box if positions are written");
847
848 if (write_data(tng, step, elapsedSeconds, reinterpret_cast<const real*>(x), 3,
849 TNG_TRAJ_POSITIONS, "POSITIONS", TNG_PARTICLE_BLOCK_DATA, compression)
850 != TNG_SUCCESS)
851 {
852 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
853 }
854 }
855
856 if (v)
857 {
858 if (write_data(tng, step, elapsedSeconds, reinterpret_cast<const real*>(v), 3,
859 TNG_TRAJ_VELOCITIES, "VELOCITIES", TNG_PARTICLE_BLOCK_DATA, compression)
860 != TNG_SUCCESS)
861 {
862 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
863 }
864 }
865
866 if (f)
867 {
868 /* TNG-MF1 compression only compresses positions and velocities. Use lossless
869 * compression for forces regardless of output mode */
870 if (write_data(tng, step, elapsedSeconds, reinterpret_cast<const real*>(f), 3,
871 TNG_TRAJ_FORCES, "FORCES", TNG_PARTICLE_BLOCK_DATA, TNG_GZIP_COMPRESSION)
872 != TNG_SUCCESS)
873 {
874 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
875 }
876 }
877
878 if (box)
879 {
880 /* TNG-MF1 compression only compresses positions and velocities. Use lossless
881 * compression for box shape regardless of output mode */
882 if (write_data(tng, step, elapsedSeconds, reinterpret_cast<const real*>(box), 9,
883 TNG_TRAJ_BOX_SHAPE, "BOX SHAPE", TNG_NON_PARTICLE_BLOCK_DATA, TNG_GZIP_COMPRESSION)
884 != TNG_SUCCESS)
885 {
886 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
887 }
888 }
889
890 if (lambda >= 0)
891 {
892 /* TNG-MF1 compression only compresses positions and velocities. Use lossless
893 * compression for lambda regardless of output mode */
894 if (write_data(tng, step, elapsedSeconds, reinterpret_cast<const real*>(&lambda), 1,
895 TNG_GMX_LAMBDA, "LAMBDAS", TNG_NON_PARTICLE_BLOCK_DATA, TNG_GZIP_COMPRESSION)
896 != TNG_SUCCESS)
897 {
898 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
899 }
900 }
901
902 // Update the data in the wrapper
903
904 gmx_tng->lastStepDataIsValid = true;
905 gmx_tng->lastStep = step;
906 gmx_tng->lastTimeDataIsValid = true;
907 gmx_tng->lastTime = elapsedSeconds;
908 #else
909 GMX_UNUSED_VALUE(gmx_tng);
910 GMX_UNUSED_VALUE(bUseLossyCompression);
911 GMX_UNUSED_VALUE(step);
912 GMX_UNUSED_VALUE(elapsedPicoSeconds);
913 GMX_UNUSED_VALUE(lambda);
914 GMX_UNUSED_VALUE(box);
915 GMX_UNUSED_VALUE(nAtoms);
916 GMX_UNUSED_VALUE(x);
917 GMX_UNUSED_VALUE(v);
918 GMX_UNUSED_VALUE(f);
919 #endif
920 }
921
fflush_tng(gmx_tng_trajectory_t gmx_tng)922 void fflush_tng(gmx_tng_trajectory_t gmx_tng)
923 {
924 #if GMX_USE_TNG
925 if (!gmx_tng)
926 {
927 return;
928 }
929 tng_frame_set_premature_write(gmx_tng->tng, TNG_USE_HASH);
930 #else
931 GMX_UNUSED_VALUE(gmx_tng);
932 #endif
933 }
934
gmx_tng_get_time_of_final_frame(gmx_tng_trajectory_t gmx_tng)935 float gmx_tng_get_time_of_final_frame(gmx_tng_trajectory_t gmx_tng)
936 {
937 #if GMX_USE_TNG
938 int64_t nFrames;
939 double time;
940 float fTime;
941 tng_trajectory_t tng = gmx_tng->tng;
942
943 tng_num_frames_get(tng, &nFrames);
944 tng_util_time_of_frame_get(tng, nFrames - 1, &time);
945
946 fTime = time / PICO;
947 return fTime;
948 #else
949 GMX_UNUSED_VALUE(gmx_tng);
950 return -1.0;
951 #endif
952 }
953
gmx_prepare_tng_writing(const char * filename,char mode,gmx_tng_trajectory_t * gmx_tng_input,gmx_tng_trajectory_t * gmx_tng_output,int nAtoms,const gmx_mtop_t * mtop,gmx::ArrayRef<const int> index,const char * indexGroupName)954 void gmx_prepare_tng_writing(const char* filename,
955 char mode,
956 gmx_tng_trajectory_t* gmx_tng_input,
957 gmx_tng_trajectory_t* gmx_tng_output,
958 int nAtoms,
959 const gmx_mtop_t* mtop,
960 gmx::ArrayRef<const int> index,
961 const char* indexGroupName)
962 {
963 #if GMX_USE_TNG
964 tng_trajectory_t* input = (gmx_tng_input && *gmx_tng_input) ? &(*gmx_tng_input)->tng : nullptr;
965 /* FIXME after 5.0: Currently only standard block types are read */
966 const int defaultNumIds = 5;
967 static int64_t fallbackIds[defaultNumIds] = { TNG_TRAJ_BOX_SHAPE, TNG_TRAJ_POSITIONS,
968 TNG_TRAJ_VELOCITIES, TNG_TRAJ_FORCES, TNG_GMX_LAMBDA };
969 static char fallbackNames[defaultNumIds][32] = { "BOX SHAPE", "POSITIONS", "VELOCITIES",
970 "FORCES", "LAMBDAS" };
971
972 typedef tng_function_status (*set_writing_interval_func_pointer)(
973 tng_trajectory_t, const int64_t, const int64_t, const int64_t, const char*, const char,
974 const char);
975 # if GMX_DOUBLE
976 set_writing_interval_func_pointer set_writing_interval = tng_util_generic_write_interval_double_set;
977 # else
978 set_writing_interval_func_pointer set_writing_interval = tng_util_generic_write_interval_set;
979 # endif
980
981 gmx_tng_open(filename, mode, gmx_tng_output);
982 tng_trajectory_t* output = &(*gmx_tng_output)->tng;
983
984 /* Do we have an input file in TNG format? If so, then there's
985 more data we can copy over, rather than having to improvise. */
986 if (gmx_tng_input && *gmx_tng_input)
987 {
988 /* Set parameters (compression, time per frame, molecule
989 * information, number of frames per frame set and writing
990 * intervals of positions, box shape and lambdas) of the
991 * output tng container based on their respective values int
992 * the input tng container */
993 double time, compression_precision;
994 int64_t n_frames_per_frame_set, interval = -1;
995
996 tng_compression_precision_get(*input, &compression_precision);
997 tng_compression_precision_set(*output, compression_precision);
998 // TODO make this configurable in a future version
999 char compression_type = TNG_TNG_COMPRESSION;
1000
1001 tng_molecule_system_copy(*input, *output);
1002
1003 tng_time_per_frame_get(*input, &time);
1004 tng_time_per_frame_set(*output, time);
1005 // Since we have copied the value from the input TNG we should not change it again
1006 (*gmx_tng_output)->timePerFrameIsSet = true;
1007
1008 tng_num_frames_per_frame_set_get(*input, &n_frames_per_frame_set);
1009 tng_num_frames_per_frame_set_set(*output, n_frames_per_frame_set);
1010
1011 for (int i = 0; i < defaultNumIds; i++)
1012 {
1013 if (tng_data_get_stride_length(*input, fallbackIds[i], -1, &interval) == TNG_SUCCESS)
1014 {
1015 switch (fallbackIds[i])
1016 {
1017 case TNG_TRAJ_POSITIONS:
1018 case TNG_TRAJ_VELOCITIES:
1019 set_writing_interval(*output, interval, 3, fallbackIds[i], fallbackNames[i],
1020 TNG_PARTICLE_BLOCK_DATA, compression_type);
1021 break;
1022 case TNG_TRAJ_FORCES:
1023 set_writing_interval(*output, interval, 3, fallbackIds[i], fallbackNames[i],
1024 TNG_PARTICLE_BLOCK_DATA, TNG_GZIP_COMPRESSION);
1025 break;
1026 case TNG_TRAJ_BOX_SHAPE:
1027 set_writing_interval(*output, interval, 9, fallbackIds[i], fallbackNames[i],
1028 TNG_NON_PARTICLE_BLOCK_DATA, TNG_GZIP_COMPRESSION);
1029 (*gmx_tng_output)->boxOutputInterval = interval;
1030 break;
1031 case TNG_GMX_LAMBDA:
1032 set_writing_interval(*output, interval, 1, fallbackIds[i], fallbackNames[i],
1033 TNG_NON_PARTICLE_BLOCK_DATA, TNG_GZIP_COMPRESSION);
1034 (*gmx_tng_output)->lambdaOutputInterval = interval;
1035 break;
1036 default: continue;
1037 }
1038 }
1039 }
1040 }
1041 else
1042 {
1043 /* TODO after trjconv is modularized: fix this so the user can
1044 change precision when they are doing an operation where
1045 this makes sense, and not otherwise.
1046
1047 char compression = bUseLossyCompression ? TNG_TNG_COMPRESSION : TNG_GZIP_COMPRESSION;
1048 gmx_tng_set_compression_precision(*output, ndec2prec(nDecimalsOfPrecision));
1049 */
1050 gmx_tng_add_mtop(*gmx_tng_output, mtop);
1051 tng_num_frames_per_frame_set_set(*output, 1);
1052 }
1053
1054 if ((!index.empty()) && nAtoms > 0)
1055 {
1056 gmx_tng_setup_atom_subgroup(*gmx_tng_output, index, indexGroupName);
1057 }
1058
1059 /* If for some reason there are more requested atoms than there are atoms in the
1060 * molecular system create a number of implicit atoms (without atom data) to
1061 * compensate for that. */
1062 if (nAtoms >= 0)
1063 {
1064 tng_implicit_num_particles_set(*output, nAtoms);
1065 }
1066 #else
1067 GMX_UNUSED_VALUE(filename);
1068 GMX_UNUSED_VALUE(mode);
1069 GMX_UNUSED_VALUE(gmx_tng_input);
1070 GMX_UNUSED_VALUE(gmx_tng_output);
1071 GMX_UNUSED_VALUE(nAtoms);
1072 GMX_UNUSED_VALUE(mtop);
1073 GMX_UNUSED_VALUE(index);
1074 GMX_UNUSED_VALUE(indexGroupName);
1075 #endif
1076 }
1077
gmx_write_tng_from_trxframe(gmx_tng_trajectory_t gmx_tng_output,const t_trxframe * frame,int natoms)1078 void gmx_write_tng_from_trxframe(gmx_tng_trajectory_t gmx_tng_output, const t_trxframe* frame, int natoms)
1079 {
1080 #if GMX_USE_TNG
1081 if (natoms < 0)
1082 {
1083 natoms = frame->natoms;
1084 }
1085 gmx_fwrite_tng(gmx_tng_output, TRUE, frame->step, frame->time, 0, frame->box, natoms, frame->x,
1086 frame->v, frame->f);
1087 #else
1088 GMX_UNUSED_VALUE(gmx_tng_output);
1089 GMX_UNUSED_VALUE(frame);
1090 GMX_UNUSED_VALUE(natoms);
1091 #endif
1092 }
1093
1094 namespace
1095 {
1096
1097 #if GMX_USE_TNG
convert_array_to_real_array(void * from,real * to,const float fact,const int nAtoms,const int nValues,const char datatype)1098 void convert_array_to_real_array(void* from,
1099 real* to,
1100 const float fact,
1101 const int nAtoms,
1102 const int nValues,
1103 const char datatype)
1104 {
1105 int i, j;
1106
1107 const bool useDouble = GMX_DOUBLE;
1108 switch (datatype)
1109 {
1110 case TNG_FLOAT_DATA:
1111 if (!useDouble)
1112 {
1113 if (fact == 1)
1114 {
1115 memcpy(to, from, nValues * sizeof(real) * nAtoms);
1116 }
1117 else
1118 {
1119 for (i = 0; i < nAtoms; i++)
1120 {
1121 for (j = 0; j < nValues; j++)
1122 {
1123 to[i * nValues + j] = reinterpret_cast<float*>(from)[i * nValues + j] * fact;
1124 }
1125 }
1126 }
1127 }
1128 else
1129 {
1130 for (i = 0; i < nAtoms; i++)
1131 {
1132 for (j = 0; j < nValues; j++)
1133 {
1134 to[i * nValues + j] = reinterpret_cast<float*>(from)[i * nValues + j] * fact;
1135 }
1136 }
1137 }
1138 break;
1139 case TNG_INT_DATA:
1140 for (i = 0; i < nAtoms; i++)
1141 {
1142 for (j = 0; j < nValues; j++)
1143 {
1144 to[i * nValues + j] = reinterpret_cast<int64_t*>(from)[i * nValues + j] * fact;
1145 }
1146 }
1147 break;
1148 case TNG_DOUBLE_DATA:
1149 if (sizeof(real) == sizeof(double))
1150 {
1151 if (fact == 1)
1152 {
1153 memcpy(to, from, nValues * sizeof(real) * nAtoms);
1154 }
1155 else
1156 {
1157 for (i = 0; i < nAtoms; i++)
1158 {
1159 for (j = 0; j < nValues; j++)
1160 {
1161 to[i * nValues + j] = reinterpret_cast<double*>(from)[i * nValues + j] * fact;
1162 }
1163 }
1164 }
1165 }
1166 else
1167 {
1168 for (i = 0; i < nAtoms; i++)
1169 {
1170 for (j = 0; j < nValues; j++)
1171 {
1172 to[i * nValues + j] = reinterpret_cast<double*>(from)[i * nValues + j] * fact;
1173 }
1174 }
1175 }
1176 break;
1177 default: gmx_incons("Illegal datatype when converting values to a real array!");
1178 }
1179 }
1180
getDistanceScaleFactor(gmx_tng_trajectory_t in)1181 real getDistanceScaleFactor(gmx_tng_trajectory_t in)
1182 {
1183 int64_t exp = -1;
1184 real distanceScaleFactor;
1185
1186 // TODO Hopefully, TNG 2.0 will do this kind of thing for us
1187 tng_distance_unit_exponential_get(in->tng, &exp);
1188
1189 // GROMACS expects distances in nm
1190 switch (exp)
1191 {
1192 case 9: distanceScaleFactor = NANO / NANO; break;
1193 case 10: distanceScaleFactor = NANO / ANGSTROM; break;
1194 default: distanceScaleFactor = pow(10.0, exp + 9.0);
1195 }
1196
1197 return distanceScaleFactor;
1198 }
1199 #endif
1200
1201 } // namespace
1202
gmx_tng_setup_atom_subgroup(gmx_tng_trajectory_t gmx_tng,gmx::ArrayRef<const int> ind,const char * name)1203 void gmx_tng_setup_atom_subgroup(gmx_tng_trajectory_t gmx_tng, gmx::ArrayRef<const int> ind, const char* name)
1204 {
1205 #if GMX_USE_TNG
1206 int64_t nAtoms, cnt, nMols;
1207 tng_molecule_t mol, iterMol;
1208 tng_chain_t chain;
1209 tng_residue_t res;
1210 tng_atom_t atom;
1211 tng_function_status stat;
1212 tng_trajectory_t tng = gmx_tng->tng;
1213
1214 tng_num_particles_get(tng, &nAtoms);
1215
1216 if (nAtoms == ind.ssize())
1217 {
1218 return;
1219 }
1220
1221 stat = tng_molecule_find(tng, name, -1, &mol);
1222 if (stat == TNG_SUCCESS)
1223 {
1224 tng_molecule_num_atoms_get(tng, mol, &nAtoms);
1225 tng_molecule_cnt_get(tng, mol, &cnt);
1226 if (nAtoms == ind.ssize())
1227 {
1228 stat = TNG_SUCCESS;
1229 }
1230 else
1231 {
1232 stat = TNG_FAILURE;
1233 }
1234 }
1235 if (stat == TNG_FAILURE)
1236 {
1237 /* The indexed atoms are added to one separate molecule. */
1238 tng_molecule_alloc(tng, &mol);
1239 tng_molecule_name_set(tng, mol, name);
1240 tng_molecule_chain_add(tng, mol, "", &chain);
1241
1242 for (gmx::index i = 0; i < ind.ssize(); i++)
1243 {
1244 char temp_name[256], temp_type[256];
1245
1246 /* Try to retrieve the residue name of the atom */
1247 stat = tng_residue_name_of_particle_nr_get(tng, ind[i], temp_name, 256);
1248 if (stat != TNG_SUCCESS)
1249 {
1250 temp_name[0] = '\0';
1251 }
1252 /* Check if the molecule of the selection already contains this residue */
1253 if (tng_chain_residue_find(tng, chain, temp_name, -1, &res) != TNG_SUCCESS)
1254 {
1255 tng_chain_residue_add(tng, chain, temp_name, &res);
1256 }
1257 /* Try to find the original name and type of the atom */
1258 stat = tng_atom_name_of_particle_nr_get(tng, ind[i], temp_name, 256);
1259 if (stat != TNG_SUCCESS)
1260 {
1261 temp_name[0] = '\0';
1262 }
1263 stat = tng_atom_type_of_particle_nr_get(tng, ind[i], temp_type, 256);
1264 if (stat != TNG_SUCCESS)
1265 {
1266 temp_type[0] = '\0';
1267 }
1268 tng_residue_atom_w_id_add(tng, res, temp_name, temp_type, ind[i], &atom);
1269 }
1270 tng_molecule_existing_add(tng, &mol);
1271 }
1272 /* Set the count of the molecule containing the selected atoms to 1 and all
1273 * other molecules to 0 */
1274 tng_molecule_cnt_set(tng, mol, 1);
1275 tng_num_molecule_types_get(tng, &nMols);
1276 for (int64_t k = 0; k < nMols; k++)
1277 {
1278 tng_molecule_of_index_get(tng, k, &iterMol);
1279 if (iterMol == mol)
1280 {
1281 continue;
1282 }
1283 tng_molecule_cnt_set(tng, iterMol, 0);
1284 }
1285 #else
1286 GMX_UNUSED_VALUE(gmx_tng);
1287 GMX_UNUSED_VALUE(ind);
1288 GMX_UNUSED_VALUE(name);
1289 #endif
1290 }
1291
1292 /* TODO: If/when TNG acquires the ability to copy data blocks without
1293 * uncompressing them, then this implemenation should be reconsidered.
1294 * Ideally, gmx trjconv -f a.tng -o b.tng -b 10 -e 20 would be fast
1295 * and lose no information. */
gmx_read_next_tng_frame(gmx_tng_trajectory_t gmx_tng_input,t_trxframe * fr,int64_t * requestedIds,int numRequestedIds)1296 gmx_bool gmx_read_next_tng_frame(gmx_tng_trajectory_t gmx_tng_input,
1297 t_trxframe* fr,
1298 int64_t* requestedIds,
1299 int numRequestedIds)
1300 {
1301 #if GMX_USE_TNG
1302 tng_trajectory_t input = gmx_tng_input->tng;
1303 gmx_bool bOK = TRUE;
1304 tng_function_status stat;
1305 int64_t numberOfAtoms = -1, frameNumber = -1;
1306 int64_t nBlocks, blockId, *blockIds = nullptr, codecId;
1307 char datatype = -1;
1308 void* values = nullptr;
1309 double frameTime = -1.0;
1310 int size, blockDependency;
1311 double prec;
1312 const int defaultNumIds = 5;
1313 static int64_t fallbackRequestedIds[defaultNumIds] = { TNG_TRAJ_BOX_SHAPE, TNG_TRAJ_POSITIONS,
1314 TNG_TRAJ_VELOCITIES, TNG_TRAJ_FORCES,
1315 TNG_GMX_LAMBDA };
1316
1317
1318 fr->bStep = FALSE;
1319 fr->bTime = FALSE;
1320 fr->bLambda = FALSE;
1321 fr->bAtoms = FALSE;
1322 fr->bPrec = FALSE;
1323 fr->bX = FALSE;
1324 fr->bV = FALSE;
1325 fr->bF = FALSE;
1326 fr->bBox = FALSE;
1327
1328 /* If no specific IDs were requested read all block types that can
1329 * currently be interpreted */
1330 if (!requestedIds || numRequestedIds == 0)
1331 {
1332 numRequestedIds = defaultNumIds;
1333 requestedIds = fallbackRequestedIds;
1334 }
1335
1336 stat = tng_num_particles_get(input, &numberOfAtoms);
1337 if (stat != TNG_SUCCESS)
1338 {
1339 gmx_file("Cannot determine number of atoms from TNG file.");
1340 }
1341 fr->natoms = numberOfAtoms;
1342
1343 bool nextFrameExists = gmx_get_tng_data_block_types_of_next_frame(
1344 gmx_tng_input, fr->step, numRequestedIds, requestedIds, &frameNumber, &nBlocks, &blockIds);
1345 gmx::unique_cptr<int64_t, gmx::free_wrapper> blockIdsGuard(blockIds);
1346 if (!nextFrameExists)
1347 {
1348 return FALSE;
1349 }
1350
1351 if (nBlocks == 0)
1352 {
1353 return FALSE;
1354 }
1355
1356 for (int64_t i = 0; i < nBlocks; i++)
1357 {
1358 blockId = blockIds[i];
1359 tng_data_block_dependency_get(input, blockId, &blockDependency);
1360 if (blockDependency & TNG_PARTICLE_DEPENDENT)
1361 {
1362 stat = tng_util_particle_data_next_frame_read(input, blockId, &values, &datatype,
1363 &frameNumber, &frameTime);
1364 }
1365 else
1366 {
1367 stat = tng_util_non_particle_data_next_frame_read(input, blockId, &values, &datatype,
1368 &frameNumber, &frameTime);
1369 }
1370 if (stat == TNG_CRITICAL)
1371 {
1372 gmx_file("Cannot read positions from TNG file.");
1373 return FALSE;
1374 }
1375 else if (stat == TNG_FAILURE)
1376 {
1377 continue;
1378 }
1379 switch (blockId)
1380 {
1381 case TNG_TRAJ_BOX_SHAPE:
1382 switch (datatype)
1383 {
1384 case TNG_INT_DATA: size = sizeof(int64_t); break;
1385 case TNG_FLOAT_DATA: size = sizeof(float); break;
1386 case TNG_DOUBLE_DATA: size = sizeof(double); break;
1387 default: gmx_incons("Illegal datatype of box shape values!");
1388 }
1389 for (int i = 0; i < DIM; i++)
1390 {
1391 convert_array_to_real_array(reinterpret_cast<char*>(values) + size * i * DIM,
1392 reinterpret_cast<real*>(fr->box[i]),
1393 getDistanceScaleFactor(gmx_tng_input), 1, DIM, datatype);
1394 }
1395 fr->bBox = TRUE;
1396 break;
1397 case TNG_TRAJ_POSITIONS:
1398 srenew(fr->x, fr->natoms);
1399 convert_array_to_real_array(values, reinterpret_cast<real*>(fr->x),
1400 getDistanceScaleFactor(gmx_tng_input), fr->natoms, DIM,
1401 datatype);
1402 fr->bX = TRUE;
1403 tng_util_frame_current_compression_get(input, blockId, &codecId, &prec);
1404 /* This must be updated if/when more lossy compression methods are added */
1405 if (codecId == TNG_TNG_COMPRESSION)
1406 {
1407 fr->prec = prec;
1408 fr->bPrec = TRUE;
1409 }
1410 break;
1411 case TNG_TRAJ_VELOCITIES:
1412 srenew(fr->v, fr->natoms);
1413 convert_array_to_real_array(values, reinterpret_cast<real*>(fr->v),
1414 getDistanceScaleFactor(gmx_tng_input), fr->natoms, DIM,
1415 datatype);
1416 fr->bV = TRUE;
1417 tng_util_frame_current_compression_get(input, blockId, &codecId, &prec);
1418 /* This must be updated if/when more lossy compression methods are added */
1419 if (codecId == TNG_TNG_COMPRESSION)
1420 {
1421 fr->prec = prec;
1422 fr->bPrec = TRUE;
1423 }
1424 break;
1425 case TNG_TRAJ_FORCES:
1426 srenew(fr->f, fr->natoms);
1427 convert_array_to_real_array(values, reinterpret_cast<real*>(fr->f),
1428 getDistanceScaleFactor(gmx_tng_input), fr->natoms, DIM,
1429 datatype);
1430 fr->bF = TRUE;
1431 break;
1432 case TNG_GMX_LAMBDA:
1433 switch (datatype)
1434 {
1435 case TNG_FLOAT_DATA: fr->lambda = *(reinterpret_cast<float*>(values)); break;
1436 case TNG_DOUBLE_DATA: fr->lambda = *(reinterpret_cast<double*>(values)); break;
1437 default: gmx_incons("Illegal datatype lambda value!");
1438 }
1439 fr->bLambda = TRUE;
1440 break;
1441 default:
1442 gmx_warning(
1443 "Illegal block type! Currently GROMACS tools can only handle certain data "
1444 "types. Skipping block.");
1445 }
1446 /* values does not have to be freed before reading next frame. It will
1447 * be reallocated if it is not NULL. */
1448 }
1449
1450 fr->step = frameNumber;
1451 fr->bStep = TRUE;
1452
1453 // Convert the time to ps
1454 fr->time = frameTime / PICO;
1455 fr->bTime = (frameTime > 0);
1456
1457 // TODO This does not leak, but is not exception safe.
1458 /* values must be freed before leaving this function */
1459 sfree(values);
1460
1461 return bOK;
1462 #else
1463 GMX_UNUSED_VALUE(gmx_tng_input);
1464 GMX_UNUSED_VALUE(fr);
1465 GMX_UNUSED_VALUE(requestedIds);
1466 GMX_UNUSED_VALUE(numRequestedIds);
1467 return FALSE;
1468 #endif
1469 }
1470
gmx_print_tng_molecule_system(gmx_tng_trajectory_t gmx_tng_input,FILE * stream)1471 void gmx_print_tng_molecule_system(gmx_tng_trajectory_t gmx_tng_input, FILE* stream)
1472 {
1473 #if GMX_USE_TNG
1474 int64_t nMolecules, nChains, nResidues, nAtoms, nFramesRead;
1475 int64_t strideLength, nParticlesRead, nValuesPerFrameRead, *molCntList;
1476 tng_molecule_t molecule;
1477 tng_chain_t chain;
1478 tng_residue_t residue;
1479 tng_atom_t atom;
1480 tng_function_status stat;
1481 char str[256];
1482 char varNAtoms;
1483 char datatype;
1484 void* data = nullptr;
1485 std::vector<real> atomCharges;
1486 std::vector<real> atomMasses;
1487 tng_trajectory_t input = gmx_tng_input->tng;
1488
1489 tng_num_molecule_types_get(input, &nMolecules);
1490 tng_molecule_cnt_list_get(input, &molCntList);
1491 /* Can the number of particles change in the trajectory or is it constant? */
1492 tng_num_particles_variable_get(input, &varNAtoms);
1493
1494 for (int64_t i = 0; i < nMolecules; i++)
1495 {
1496 tng_molecule_of_index_get(input, i, &molecule);
1497 tng_molecule_name_get(input, molecule, str, 256);
1498 if (varNAtoms == TNG_CONSTANT_N_ATOMS)
1499 {
1500 if (static_cast<int>(molCntList[i]) == 0)
1501 {
1502 continue;
1503 }
1504 fprintf(stream, "Molecule: %s, count: %d\n", str, static_cast<int>(molCntList[i]));
1505 }
1506 else
1507 {
1508 fprintf(stream, "Molecule: %s\n", str);
1509 }
1510 tng_molecule_num_chains_get(input, molecule, &nChains);
1511 if (nChains > 0)
1512 {
1513 for (int64_t j = 0; j < nChains; j++)
1514 {
1515 tng_molecule_chain_of_index_get(input, molecule, j, &chain);
1516 tng_chain_name_get(input, chain, str, 256);
1517 fprintf(stream, "\tChain: %s\n", str);
1518 tng_chain_num_residues_get(input, chain, &nResidues);
1519 for (int64_t k = 0; k < nResidues; k++)
1520 {
1521 tng_chain_residue_of_index_get(input, chain, k, &residue);
1522 tng_residue_name_get(input, residue, str, 256);
1523 fprintf(stream, "\t\tResidue: %s\n", str);
1524 tng_residue_num_atoms_get(input, residue, &nAtoms);
1525 for (int64_t l = 0; l < nAtoms; l++)
1526 {
1527 tng_residue_atom_of_index_get(input, residue, l, &atom);
1528 tng_atom_name_get(input, atom, str, 256);
1529 fprintf(stream, "\t\t\tAtom: %s", str);
1530 tng_atom_type_get(input, atom, str, 256);
1531 fprintf(stream, " (%s)\n", str);
1532 }
1533 }
1534 }
1535 }
1536 /* It is possible to have a molecule without chains, in which case
1537 * residues in the molecule can be iterated through without going
1538 * through chains. */
1539 else
1540 {
1541 tng_molecule_num_residues_get(input, molecule, &nResidues);
1542 if (nResidues > 0)
1543 {
1544 for (int64_t k = 0; k < nResidues; k++)
1545 {
1546 tng_molecule_residue_of_index_get(input, molecule, k, &residue);
1547 tng_residue_name_get(input, residue, str, 256);
1548 fprintf(stream, "\t\tResidue: %s\n", str);
1549 tng_residue_num_atoms_get(input, residue, &nAtoms);
1550 for (int64_t l = 0; l < nAtoms; l++)
1551 {
1552 tng_residue_atom_of_index_get(input, residue, l, &atom);
1553 tng_atom_name_get(input, atom, str, 256);
1554 fprintf(stream, "\t\t\tAtom: %s", str);
1555 tng_atom_type_get(input, atom, str, 256);
1556 fprintf(stream, " (%s)\n", str);
1557 }
1558 }
1559 }
1560 else
1561 {
1562 tng_molecule_num_atoms_get(input, molecule, &nAtoms);
1563 for (int64_t l = 0; l < nAtoms; l++)
1564 {
1565 tng_molecule_atom_of_index_get(input, molecule, l, &atom);
1566 tng_atom_name_get(input, atom, str, 256);
1567 fprintf(stream, "\t\t\tAtom: %s", str);
1568 tng_atom_type_get(input, atom, str, 256);
1569 fprintf(stream, " (%s)\n", str);
1570 }
1571 }
1572 }
1573 }
1574
1575 tng_num_particles_get(input, &nAtoms);
1576 stat = tng_particle_data_vector_get(input, TNG_TRAJ_PARTIAL_CHARGES, &data, &nFramesRead,
1577 &strideLength, &nParticlesRead, &nValuesPerFrameRead, &datatype);
1578 if (stat == TNG_SUCCESS)
1579 {
1580 atomCharges.resize(nAtoms);
1581 convert_array_to_real_array(data, atomCharges.data(), 1, nAtoms, 1, datatype);
1582
1583 fprintf(stream, "Atom Charges (%d):\n", int(nAtoms));
1584 for (int64_t i = 0; i < nAtoms; i += 10)
1585 {
1586 fprintf(stream, "Atom Charges [%8d-]=[", int(i));
1587 for (int64_t j = 0; (j < 10 && i + j < nAtoms); j++)
1588 {
1589 fprintf(stream, " %12.5e", atomCharges[i + j]);
1590 }
1591 fprintf(stream, "]\n");
1592 }
1593 }
1594
1595 stat = tng_particle_data_vector_get(input, TNG_TRAJ_MASSES, &data, &nFramesRead, &strideLength,
1596 &nParticlesRead, &nValuesPerFrameRead, &datatype);
1597 if (stat == TNG_SUCCESS)
1598 {
1599 atomMasses.resize(nAtoms);
1600 convert_array_to_real_array(data, atomMasses.data(), 1, nAtoms, 1, datatype);
1601
1602 fprintf(stream, "Atom Masses (%d):\n", int(nAtoms));
1603 for (int64_t i = 0; i < nAtoms; i += 10)
1604 {
1605 fprintf(stream, "Atom Masses [%8d-]=[", int(i));
1606 for (int64_t j = 0; (j < 10 && i + j < nAtoms); j++)
1607 {
1608 fprintf(stream, " %12.5e", atomMasses[i + j]);
1609 }
1610 fprintf(stream, "]\n");
1611 }
1612 }
1613
1614 sfree(data);
1615 #else
1616 GMX_UNUSED_VALUE(gmx_tng_input);
1617 GMX_UNUSED_VALUE(stream);
1618 #endif
1619 }
1620
gmx_get_tng_data_block_types_of_next_frame(gmx_tng_trajectory_t gmx_tng_input,int frame,int nRequestedIds,int64_t * requestedIds,int64_t * nextFrame,int64_t * nBlocks,int64_t ** blockIds)1621 gmx_bool gmx_get_tng_data_block_types_of_next_frame(gmx_tng_trajectory_t gmx_tng_input,
1622 int frame,
1623 int nRequestedIds,
1624 int64_t* requestedIds,
1625 int64_t* nextFrame,
1626 int64_t* nBlocks,
1627 int64_t** blockIds)
1628 {
1629 #if GMX_USE_TNG
1630 tng_function_status stat;
1631 tng_trajectory_t input = gmx_tng_input->tng;
1632
1633 stat = tng_util_trajectory_next_frame_present_data_blocks_find(
1634 input, frame, nRequestedIds, requestedIds, nextFrame, nBlocks, blockIds);
1635
1636 if (stat == TNG_CRITICAL)
1637 {
1638 gmx_file("Cannot read TNG file. Cannot find data blocks of next frame.");
1639 }
1640 else if (stat == TNG_FAILURE)
1641 {
1642 return FALSE;
1643 }
1644 return TRUE;
1645 #else
1646 GMX_UNUSED_VALUE(gmx_tng_input);
1647 GMX_UNUSED_VALUE(frame);
1648 GMX_UNUSED_VALUE(nRequestedIds);
1649 GMX_UNUSED_VALUE(requestedIds);
1650 GMX_UNUSED_VALUE(nextFrame);
1651 GMX_UNUSED_VALUE(nBlocks);
1652 GMX_UNUSED_VALUE(blockIds);
1653 return FALSE;
1654 #endif
1655 }
1656
gmx_get_tng_data_next_frame_of_block_type(gmx_tng_trajectory_t gmx_tng_input,int64_t blockId,real ** values,int64_t * frameNumber,double * frameTime,int64_t * nValuesPerFrame,int64_t * nAtoms,real * prec,char * name,int maxLen,gmx_bool * bOK)1657 gmx_bool gmx_get_tng_data_next_frame_of_block_type(gmx_tng_trajectory_t gmx_tng_input,
1658 int64_t blockId,
1659 real** values,
1660 int64_t* frameNumber,
1661 double* frameTime,
1662 int64_t* nValuesPerFrame,
1663 int64_t* nAtoms,
1664 real* prec,
1665 char* name,
1666 int maxLen,
1667 gmx_bool* bOK)
1668 {
1669 #if GMX_USE_TNG
1670 tng_function_status stat;
1671 char datatype = -1;
1672 int64_t codecId;
1673 int blockDependency;
1674 void* data = nullptr;
1675 double localPrec;
1676 tng_trajectory_t input = gmx_tng_input->tng;
1677
1678 stat = tng_data_block_name_get(input, blockId, name, maxLen);
1679 if (stat != TNG_SUCCESS)
1680 {
1681 gmx_file("Cannot read next frame of TNG file");
1682 }
1683 stat = tng_data_block_dependency_get(input, blockId, &blockDependency);
1684 if (stat != TNG_SUCCESS)
1685 {
1686 gmx_file("Cannot read next frame of TNG file");
1687 }
1688 if (blockDependency & TNG_PARTICLE_DEPENDENT)
1689 {
1690 tng_num_particles_get(input, nAtoms);
1691 stat = tng_util_particle_data_next_frame_read(input, blockId, &data, &datatype, frameNumber,
1692 frameTime);
1693 }
1694 else
1695 {
1696 *nAtoms = 1; /* There are not actually any atoms, but it is used for
1697 allocating memory */
1698 stat = tng_util_non_particle_data_next_frame_read(input, blockId, &data, &datatype,
1699 frameNumber, frameTime);
1700 }
1701 if (stat == TNG_CRITICAL)
1702 {
1703 gmx_file("Cannot read next frame of TNG file");
1704 }
1705 if (stat == TNG_FAILURE)
1706 {
1707 *bOK = TRUE;
1708 return FALSE;
1709 }
1710
1711 stat = tng_data_block_num_values_per_frame_get(input, blockId, nValuesPerFrame);
1712 if (stat != TNG_SUCCESS)
1713 {
1714 gmx_file("Cannot read next frame of TNG file");
1715 }
1716 srenew(*values, sizeof(real) * *nValuesPerFrame * *nAtoms);
1717 convert_array_to_real_array(data, *values, getDistanceScaleFactor(gmx_tng_input), *nAtoms,
1718 *nValuesPerFrame, datatype);
1719
1720 tng_util_frame_current_compression_get(input, blockId, &codecId, &localPrec);
1721
1722 /* This must be updated if/when more lossy compression methods are added */
1723 if (codecId != TNG_TNG_COMPRESSION)
1724 {
1725 *prec = -1.0;
1726 }
1727 else
1728 {
1729 *prec = localPrec;
1730 }
1731
1732 sfree(data);
1733
1734 *bOK = TRUE;
1735 return TRUE;
1736 #else
1737 GMX_UNUSED_VALUE(gmx_tng_input);
1738 GMX_UNUSED_VALUE(blockId);
1739 GMX_UNUSED_VALUE(values);
1740 GMX_UNUSED_VALUE(frameNumber);
1741 GMX_UNUSED_VALUE(frameTime);
1742 GMX_UNUSED_VALUE(nValuesPerFrame);
1743 GMX_UNUSED_VALUE(nAtoms);
1744 GMX_UNUSED_VALUE(prec);
1745 GMX_UNUSED_VALUE(name);
1746 GMX_UNUSED_VALUE(maxLen);
1747 GMX_UNUSED_VALUE(bOK);
1748 return FALSE;
1749 #endif
1750 }
1751
gmx_tng_get_box_output_interval(gmx_tng_trajectory_t gmx_tng)1752 int gmx_tng_get_box_output_interval(gmx_tng_trajectory_t gmx_tng)
1753 {
1754 #if GMX_USE_TNG
1755 return gmx_tng->boxOutputInterval;
1756 #else
1757 GMX_UNUSED_VALUE(gmx_tng);
1758 return -1;
1759 #endif
1760 }
1761
gmx_tng_get_lambda_output_interval(gmx_tng_trajectory_t gmx_tng)1762 int gmx_tng_get_lambda_output_interval(gmx_tng_trajectory_t gmx_tng)
1763 {
1764 #if GMX_USE_TNG
1765 return gmx_tng->lambdaOutputInterval;
1766 #else
1767 GMX_UNUSED_VALUE(gmx_tng);
1768 return -1;
1769 #endif
1770 }
1771