2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014,2015,2018, 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.
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.
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.
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.
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.
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.
49 #include "tng/tng_io.h"
52 #include "gromacs/math/units.h"
53 #include "gromacs/math/utilities.h"
54 #include "gromacs/mdtypes/inputrec.h"
55 #include "gromacs/topology/ifunc.h"
56 #include "gromacs/topology/topology.h"
57 #include "gromacs/trajectory/trajectoryframe.h"
58 #include "gromacs/utility/basedefinitions.h"
59 #include "gromacs/utility/baseversion.h"
60 #include "gromacs/utility/fatalerror.h"
61 #include "gromacs/utility/futil.h"
62 #include "gromacs/utility/gmxassert.h"
63 #include "gromacs/utility/programcontext.h"
64 #include "gromacs/utility/smalloc.h"
65 #include "gromacs/utility/sysinfo.h"
66 #include "gromacs/utility/unique_cptr.h"
68 /*! \brief Gromacs Wrapper around tng datatype
70 * This could in principle hold any GROMACS-specific requirements not yet
71 * implemented in or not relevant to the TNG library itself. However, for now
72 * we only use it to handle some shortcomings we have discovered, where the TNG
73 * API itself is a bit fragile and can end up overwriting data if called several
74 * times with the same frame number.
75 * The logic to determine the time per step was also a bit fragile. This is not
76 * critical, but since we anyway need a wrapper for ensuring unique frame
77 * numbers, we can also use it to store the time of the first step and use that
78 * to derive a slightly better/safer estimate of the time per step.
80 * At some future point where we have a second-generation TNG API we should
81 * consider removing this again.
83 struct gmx_tng_trajectory
85 tng_trajectory_t tng
; //!< Actual TNG handle (pointer)
86 bool lastStepDataIsValid
; //!< True if lastStep has been set
87 std::int64_t lastStep
; //!< Index/step used for last frame
88 bool lastTimeDataIsValid
; //!< True if lastTime has been set
89 double lastTime
; //!< Time of last frame (TNG unit is seconds)
90 bool timePerFrameIsSet
; //!< True if we have set the time per frame
93 static const char *modeToVerb(char mode
)
108 gmx_fatal(FARGS
, "Invalid file opening mode %c", mode
);
115 void gmx_tng_open(const char *filename
,
117 gmx_tng_trajectory_t
*gmx_tng
)
120 /* First check whether we have to make a backup,
121 * only for writing, not for read or append.
125 make_backup(filename
);
128 *gmx_tng
= new gmx_tng_trajectory
;
129 (*gmx_tng
)->lastStepDataIsValid
= false;
130 (*gmx_tng
)->lastTimeDataIsValid
= false;
131 (*gmx_tng
)->timePerFrameIsSet
= false;
132 tng_trajectory_t
* tng
= &(*gmx_tng
)->tng
;
134 /* tng must not be pointing at already allocated memory.
135 * Memory will be allocated by tng_util_trajectory_open() and must
136 * later on be freed by tng_util_trajectory_close(). */
137 if (TNG_SUCCESS
!= tng_util_trajectory_open(filename
, mode
, tng
))
139 /* TNG does return more than one degree of error, but there is
140 no use case for GROMACS handling the non-fatal errors
143 "File I/O error while opening %s for %s",
148 if (mode
== 'w' || mode
== 'a')
151 gmx_gethostname(hostname
, 256);
154 tng_first_computer_name_set(*tng
, hostname
);
158 tng_last_computer_name_set(*tng
, hostname
);
161 char programInfo
[256];
162 const char *precisionString
= "";
164 precisionString
= " (double precision)";
166 sprintf(programInfo
, "%.100s %.128s%.24s",
167 gmx::getProgramContext().displayName(),
168 gmx_version(), precisionString
);
171 tng_first_program_name_set(*tng
, programInfo
);
175 tng_last_program_name_set(*tng
, programInfo
);
179 if (!gmx_getusername(username
, 256))
183 tng_first_user_name_set(*tng
, username
);
187 tng_last_user_name_set(*tng
, username
);
188 tng_file_headers_write(*tng
, TNG_USE_HASH
);
193 gmx_file("GROMACS was compiled without TNG support, cannot handle this file type");
194 GMX_UNUSED_VALUE(filename
);
195 GMX_UNUSED_VALUE(mode
);
196 GMX_UNUSED_VALUE(gmx_tng
);
200 void gmx_tng_close(gmx_tng_trajectory_t
*gmx_tng
)
202 /* We have to check that tng is set because
203 * tng_util_trajectory_close wants to return a NULL in it, and
204 * gives a fatal error if it is NULL. */
206 if (gmx_tng
== nullptr || *gmx_tng
== nullptr)
210 tng_trajectory_t
* tng
= &(*gmx_tng
)->tng
;
214 tng_util_trajectory_close(tng
);
220 GMX_UNUSED_VALUE(gmx_tng
);
225 static void addTngMoleculeFromTopology(gmx_tng_trajectory_t gmx_tng
,
226 const char *moleculeName
,
227 const t_atoms
*atoms
,
228 gmx_int64_t numMolecules
,
229 tng_molecule_t
*tngMol
)
231 tng_trajectory_t tng
= gmx_tng
->tng
;
232 tng_chain_t tngChain
= nullptr;
233 tng_residue_t tngRes
= nullptr;
235 if (tng_molecule_add(tng
, moleculeName
, tngMol
) != TNG_SUCCESS
)
237 gmx_file("Cannot add molecule to TNG molecular system.");
240 for (int atomIndex
= 0; atomIndex
< atoms
->nr
; atomIndex
++)
242 const t_atom
*at
= &atoms
->atom
[atomIndex
];
243 /* FIXME: Currently the TNG API can only add atoms belonging to a
244 * residue and chain. Wait for TNG 2.0*/
247 const t_resinfo
*resInfo
= &atoms
->resinfo
[at
->resind
];
248 char chainName
[2] = {resInfo
->chainid
, 0};
249 tng_atom_t tngAtom
= nullptr;
254 prevAtom
= &atoms
->atom
[atomIndex
- 1];
261 /* If this is the first atom or if the residue changed add the
262 * residue to the TNG molecular system. */
263 if (!prevAtom
|| resInfo
!= &atoms
->resinfo
[prevAtom
->resind
])
265 /* If this is the first atom or if the chain changed add
266 * the chain to the TNG molecular system. */
267 if (!prevAtom
|| resInfo
->chainid
!=
268 atoms
->resinfo
[prevAtom
->resind
].chainid
)
270 tng_molecule_chain_add(tng
, *tngMol
, chainName
,
273 /* FIXME: When TNG supports both residue index and residue
274 * number the latter should be used. Wait for TNG 2.0*/
275 tng_chain_residue_add(tng
, tngChain
, *resInfo
->name
, &tngRes
);
277 tng_residue_atom_add(tng
, tngRes
, *(atoms
->atomname
[atomIndex
]), *(atoms
->atomtype
[atomIndex
]), &tngAtom
);
280 tng_molecule_cnt_set(tng
, *tngMol
, numMolecules
);
283 void gmx_tng_add_mtop(gmx_tng_trajectory_t gmx_tng
,
284 const gmx_mtop_t
*mtop
)
288 std::vector
<real
> atomCharges
;
289 std::vector
<real
> atomMasses
;
290 const t_ilist
*ilist
;
294 tng_trajectory_t tng
= gmx_tng
->tng
;
298 /* No topology information available to add. */
303 datatype
= TNG_DOUBLE_DATA
;
305 datatype
= TNG_FLOAT_DATA
;
308 atomCharges
.reserve(mtop
->natoms
);
309 atomMasses
.reserve(mtop
->natoms
);
311 for (const gmx_molblock_t
&molBlock
: mtop
->molblock
)
313 tng_molecule_t tngMol
= nullptr;
314 const gmx_moltype_t
*molType
= &mtop
->moltype
[molBlock
.type
];
316 /* Add a molecule to the TNG trajectory with the same name as the
317 * current molecule. */
318 addTngMoleculeFromTopology(gmx_tng
,
324 /* Bonds have to be deduced from interactions (constraints etc). Different
325 * interactions have different sets of parameters. */
326 /* Constraints are specified using two atoms */
327 for (i
= 0; i
< F_NRE
; i
++)
331 ilist
= &molType
->ilist
[i
];
335 while (j
< ilist
->nr
)
337 tng_molecule_bond_add(tng
, tngMol
, ilist
->iatoms
[j
], ilist
->iatoms
[j
+1], &tngBond
);
343 /* Settle is described using three atoms */
344 ilist
= &molType
->ilist
[F_SETTLE
];
348 while (j
< ilist
->nr
)
350 tng_molecule_bond_add(tng
, tngMol
, ilist
->iatoms
[j
], ilist
->iatoms
[j
+1], &tngBond
);
351 tng_molecule_bond_add(tng
, tngMol
, ilist
->iatoms
[j
], ilist
->iatoms
[j
+2], &tngBond
);
355 /* First copy atom charges and masses, first atom by atom and then copy the memory for the molecule instances.
356 * FIXME: Atom B state data should also be written to TNG (v 2.0?) */
357 for (int atomCounter
= 0; atomCounter
< molType
->atoms
.nr
; atomCounter
++)
359 atomCharges
.push_back(molType
->atoms
.atom
[atomCounter
].q
);
360 atomMasses
.push_back(molType
->atoms
.atom
[atomCounter
].m
);
362 for (int molCounter
= 1; molCounter
< molBlock
.nmol
; molCounter
++)
364 std::copy_n(atomCharges
.end() - molType
->atoms
.nr
, molType
->atoms
.nr
, std::back_inserter(atomCharges
));
365 std::copy_n(atomMasses
.end() - molType
->atoms
.nr
, molType
->atoms
.nr
, std::back_inserter(atomMasses
));
368 /* Write the TNG data blocks. */
369 tng_particle_data_block_add(tng
, TNG_TRAJ_PARTIAL_CHARGES
, "PARTIAL CHARGES",
370 datatype
, TNG_NON_TRAJECTORY_BLOCK
,
371 1, 1, 1, 0, mtop
->natoms
,
372 TNG_GZIP_COMPRESSION
, atomCharges
.data());
373 tng_particle_data_block_add(tng
, TNG_TRAJ_MASSES
, "ATOM MASSES",
374 datatype
, TNG_NON_TRAJECTORY_BLOCK
,
375 1, 1, 1, 0, mtop
->natoms
,
376 TNG_GZIP_COMPRESSION
, atomMasses
.data());
379 /*! \libinternal \brief Compute greatest common divisor of n1 and n2
380 * if they are positive.
382 * If only one of n1 and n2 is positive, then return it.
383 * If neither n1 or n2 is positive, then return -1. */
385 greatest_common_divisor_if_positive(int n1
, int n2
)
389 return (0 >= n2
) ? -1 : n2
;
396 /* We have a non-trivial greatest common divisor to compute. */
397 return gmx_greatest_common_divisor(n1
, n2
);
400 /* By default try to write 100 frames (of actual output) in each frame set.
401 * This number is the number of outputs of the most frequently written data
402 * type per frame set.
403 * TODO for 5.1: Verify that 100 frames per frame set is efficient for most
404 * setups regarding compression efficiency and compression time. Make this
405 * a hidden command-line option? */
406 const int defaultFramesPerFrameSet
= 100;
408 /*! \libinternal \brief Set the number of frames per frame
409 * set according to output intervals.
410 * The default is that 100 frames are written of the data
411 * that is written most often. */
412 static void tng_set_frames_per_frame_set(gmx_tng_trajectory_t gmx_tng
,
413 const gmx_bool bUseLossyCompression
,
414 const t_inputrec
*ir
)
417 tng_trajectory_t tng
= gmx_tng
->tng
;
419 /* Set the number of frames per frame set to contain at least
420 * defaultFramesPerFrameSet of the lowest common denominator of
421 * the writing interval of positions and velocities. */
422 /* FIXME after 5.0: consider nstenergy also? */
423 if (bUseLossyCompression
)
425 gcd
= ir
->nstxout_compressed
;
429 gcd
= greatest_common_divisor_if_positive(ir
->nstxout
, ir
->nstvout
);
430 gcd
= greatest_common_divisor_if_positive(gcd
, ir
->nstfout
);
437 tng_num_frames_per_frame_set_set(tng
, gcd
* defaultFramesPerFrameSet
);
440 /*! \libinternal \brief Set the data-writing intervals, and number of
441 * frames per frame set */
442 static void set_writing_intervals(gmx_tng_trajectory_t gmx_tng
,
443 const gmx_bool bUseLossyCompression
,
444 const t_inputrec
*ir
)
446 tng_trajectory_t tng
= gmx_tng
->tng
;
448 /* Define pointers to specific writing functions depending on if we
449 * write float or double data */
450 typedef tng_function_status (*set_writing_interval_func_pointer
)(tng_trajectory_t
,
458 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_double_set
;
460 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_set
;
462 int xout
, vout
, fout
;
463 int gcd
= -1, lowest
= -1;
466 tng_set_frames_per_frame_set(gmx_tng
, bUseLossyCompression
, ir
);
468 if (bUseLossyCompression
)
470 xout
= ir
->nstxout_compressed
;
472 /* If there is no uncompressed coordinate output write forces
473 and velocities to the compressed tng file. */
484 compression
= TNG_TNG_COMPRESSION
;
491 compression
= TNG_GZIP_COMPRESSION
;
495 set_writing_interval(tng
, xout
, 3, TNG_TRAJ_POSITIONS
,
496 "POSITIONS", TNG_PARTICLE_BLOCK_DATA
,
498 /* TODO: if/when we write energies to TNG also, reconsider how
499 * and when box information is written, because GROMACS
500 * behaviour pre-5.0 was to write the box with every
501 * trajectory frame and every energy frame, and probably
502 * people depend on this. */
504 gcd
= greatest_common_divisor_if_positive(gcd
, xout
);
505 if (lowest
< 0 || xout
< lowest
)
512 set_writing_interval(tng
, vout
, 3, TNG_TRAJ_VELOCITIES
,
513 "VELOCITIES", TNG_PARTICLE_BLOCK_DATA
,
516 gcd
= greatest_common_divisor_if_positive(gcd
, vout
);
517 if (lowest
< 0 || vout
< lowest
)
524 set_writing_interval(tng
, fout
, 3, TNG_TRAJ_FORCES
,
525 "FORCES", TNG_PARTICLE_BLOCK_DATA
,
526 TNG_GZIP_COMPRESSION
);
528 gcd
= greatest_common_divisor_if_positive(gcd
, fout
);
529 if (lowest
< 0 || fout
< lowest
)
536 /* Lambdas and box shape written at an interval of the lowest common
537 denominator of other output */
538 set_writing_interval(tng
, gcd
, 1, TNG_GMX_LAMBDA
,
539 "LAMBDAS", TNG_NON_PARTICLE_BLOCK_DATA
,
540 TNG_GZIP_COMPRESSION
);
542 set_writing_interval(tng
, gcd
, 9, TNG_TRAJ_BOX_SHAPE
,
543 "BOX SHAPE", TNG_NON_PARTICLE_BLOCK_DATA
,
544 TNG_GZIP_COMPRESSION
);
545 if (gcd
< lowest
/ 10)
547 gmx_warning("The lowest common denominator of trajectory output is "
548 "every %d step(s), whereas the shortest output interval "
549 "is every %d steps.", gcd
, lowest
);
555 void gmx_tng_prepare_md_writing(gmx_tng_trajectory_t gmx_tng
,
556 const gmx_mtop_t
*mtop
,
557 const t_inputrec
*ir
)
560 gmx_tng_add_mtop(gmx_tng
, mtop
);
561 set_writing_intervals(gmx_tng
, FALSE
, ir
);
562 tng_time_per_frame_set(gmx_tng
->tng
, ir
->delta_t
* PICO
);
563 gmx_tng
->timePerFrameIsSet
= true;
565 GMX_UNUSED_VALUE(gmx_tng
);
566 GMX_UNUSED_VALUE(mtop
);
567 GMX_UNUSED_VALUE(ir
);
572 /* Check if all atoms in the molecule system are selected
573 * by a selection group of type specified by gtype. */
574 static gmx_bool
all_atoms_selected(const gmx_mtop_t
*mtop
,
577 /* Iterate through all atoms in the molecule system and
578 * check if they belong to a selection group of the
581 for (const gmx_molblock_t
&molBlock
: mtop
->molblock
)
583 const gmx_moltype_t
&molType
= mtop
->moltype
[molBlock
.type
];
584 const t_atoms
&atoms
= molType
.atoms
;
586 for (int j
= 0; j
< molBlock
.nmol
; j
++)
588 for (int atomIndex
= 0; atomIndex
< atoms
.nr
; atomIndex
++, i
++)
590 if (ggrpnr(&mtop
->groups
, gtype
, i
) != 0)
600 /* Create TNG molecules which will represent each of the selection
601 * groups for writing. But do that only if there is actually a
602 * specified selection group and it is not the whole system.
603 * TODO: Currently the only selection that is taken into account
604 * is egcCompressedX, but other selections should be added when
605 * e.g. writing energies is implemented.
607 static void add_selection_groups(gmx_tng_trajectory_t gmx_tng
,
608 const gmx_mtop_t
*mtop
)
610 const t_atoms
*atoms
;
612 const t_resinfo
*resInfo
;
613 const t_ilist
*ilist
;
616 tng_molecule_t mol
, iterMol
;
623 tng_trajectory_t tng
= gmx_tng
->tng
;
625 /* TODO: When the TNG molecules block is more flexible TNG selection
626 * groups should not need all atoms specified. It should be possible
627 * just to specify what molecules are selected (and/or which atoms in
628 * the molecule) and how many (if applicable). */
630 /* If no atoms are selected we do not need to create a
631 * TNG selection group molecule. */
632 if (mtop
->groups
.ngrpnr
[egcCompressedX
] == 0)
637 /* If all atoms are selected we do not have to create a selection
638 * group molecule in the TNG molecule system. */
639 if (all_atoms_selected(mtop
, egcCompressedX
))
644 /* The name of the TNG molecule containing the selection group is the
645 * same as the name of the selection group. */
646 nameIndex
= *mtop
->groups
.grps
[egcCompressedX
].nm_ind
;
647 groupName
= *mtop
->groups
.grpname
[nameIndex
];
649 tng_molecule_alloc(tng
, &mol
);
650 tng_molecule_name_set(tng
, mol
, groupName
);
651 tng_molecule_chain_add(tng
, mol
, "", &chain
);
653 for (const gmx_molblock_t
&molBlock
: mtop
->molblock
)
655 const gmx_moltype_t
&molType
= mtop
->moltype
[molBlock
.type
];
657 atoms
= &molType
.atoms
;
659 for (int j
= 0; j
< molBlock
.nmol
; j
++)
661 bool bAtomsAdded
= FALSE
;
662 for (int atomIndex
= 0; atomIndex
< atoms
->nr
; atomIndex
++, i
++)
667 if (ggrpnr(&mtop
->groups
, egcCompressedX
, i
) != 0)
671 at
= &atoms
->atom
[atomIndex
];
674 resInfo
= &atoms
->resinfo
[at
->resind
];
675 /* FIXME: When TNG supports both residue index and residue
676 * number the latter should be used. */
677 res_name
= *resInfo
->name
;
678 res_id
= at
->resind
+ 1;
682 res_name
= (char *)"";
685 if (tng_chain_residue_find(tng
, chain
, res_name
, res_id
, &res
)
688 /* Since there is ONE chain for selection groups do not keep the
689 * original residue IDs - otherwise there might be conflicts. */
690 tng_chain_residue_add(tng
, chain
, res_name
, &res
);
692 tng_residue_atom_w_id_add(tng
, res
, *(atoms
->atomname
[atomIndex
]),
693 *(atoms
->atomtype
[atomIndex
]),
694 atom_offset
+ atomIndex
, &atom
);
700 for (int k
= 0; k
< F_NRE
; k
++)
704 ilist
= &molType
.ilist
[k
];
708 while (l
< ilist
->nr
)
711 atom1
= ilist
->iatoms
[l
] + atom_offset
;
712 atom2
= ilist
->iatoms
[l
+1] + atom_offset
;
713 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom1
) == 0 &&
714 ggrpnr(&mtop
->groups
, egcCompressedX
, atom2
) == 0)
716 tng_molecule_bond_add(tng
, mol
, ilist
->iatoms
[l
],
717 ilist
->iatoms
[l
+1], &tngBond
);
724 /* Settle is described using three atoms */
725 ilist
= &molType
.ilist
[F_SETTLE
];
729 while (l
< ilist
->nr
)
731 int atom1
, atom2
, atom3
;
732 atom1
= ilist
->iatoms
[l
] + atom_offset
;
733 atom2
= ilist
->iatoms
[l
+1] + atom_offset
;
734 atom3
= ilist
->iatoms
[l
+2] + atom_offset
;
735 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom1
) == 0)
737 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom2
) == 0)
739 tng_molecule_bond_add(tng
, mol
, atom1
,
742 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom3
) == 0)
744 tng_molecule_bond_add(tng
, mol
, atom1
,
752 atom_offset
+= atoms
->nr
;
755 tng_molecule_existing_add(tng
, &mol
);
756 tng_molecule_cnt_set(tng
, mol
, 1);
757 tng_num_molecule_types_get(tng
, &nMols
);
758 for (gmx_int64_t k
= 0; k
< nMols
; k
++)
760 tng_molecule_of_index_get(tng
, k
, &iterMol
);
765 tng_molecule_cnt_set(tng
, iterMol
, 0);
770 void gmx_tng_set_compression_precision(gmx_tng_trajectory_t gmx_tng
,
774 tng_compression_precision_set(gmx_tng
->tng
, prec
);
776 GMX_UNUSED_VALUE(gmx_tng
);
777 GMX_UNUSED_VALUE(prec
);
781 void gmx_tng_prepare_low_prec_writing(gmx_tng_trajectory_t gmx_tng
,
782 const gmx_mtop_t
*mtop
,
783 const t_inputrec
*ir
)
786 gmx_tng_add_mtop(gmx_tng
, mtop
);
787 add_selection_groups(gmx_tng
, mtop
);
788 set_writing_intervals(gmx_tng
, TRUE
, ir
);
789 tng_time_per_frame_set(gmx_tng
->tng
, ir
->delta_t
* PICO
);
790 gmx_tng
->timePerFrameIsSet
= true;
791 gmx_tng_set_compression_precision(gmx_tng
, ir
->x_compression_precision
);
793 GMX_UNUSED_VALUE(gmx_tng
);
794 GMX_UNUSED_VALUE(mtop
);
795 GMX_UNUSED_VALUE(ir
);
799 void gmx_fwrite_tng(gmx_tng_trajectory_t gmx_tng
,
800 const gmx_bool bUseLossyCompression
,
802 real elapsedPicoSeconds
,
811 typedef tng_function_status (*write_data_func_pointer
)(tng_trajectory_t
,
821 static write_data_func_pointer write_data
= tng_util_generic_with_time_double_write
;
823 static write_data_func_pointer write_data
= tng_util_generic_with_time_write
;
825 double elapsedSeconds
= elapsedPicoSeconds
* PICO
;
826 gmx_int64_t nParticles
;
832 /* This function might get called when the type of the
833 compressed trajectory is actually XTC. So we exit and move
837 tng_trajectory_t tng
= gmx_tng
->tng
;
839 // While the GROMACS interface to this routine specifies 'step', TNG itself
840 // only uses 'frame index' internally, although it suggests that it's a good
841 // idea to let this represent the step rather than just frame numbers.
843 // However, the way the GROMACS interface works, there are cases where
844 // the step is simply not set, so all frames rather have step=0.
845 // If we call the lower-level TNG routines with the same frame index
846 // (which is set from the step) multiple times, things will go very
847 // wrong (overwritten data).
849 // To avoid this, we require the step value to always be larger than the
850 // previous value, and if this is not true we simply set it to a value
851 // one unit larger than the previous step.
853 // Note: We do allow the user to specify any crazy value the want for the
854 // first step. If it is "not set", GROMACS will have used the value 0.
855 if (gmx_tng
->lastStepDataIsValid
&& step
<= gmx_tng
->lastStep
)
857 step
= gmx_tng
->lastStep
+ 1;
860 // Now that we have fixed the potentially incorrect step, we can also set
861 // the time per frame if it was not already done, and we have
862 // step/time information from the last frame so it is possible to calculate it.
863 // Note that we do allow the time to be the same for all steps, or even
864 // decreasing. In that case we will simply say the time per step is 0
865 // or negative. A bit strange, but a correct representation of the data :-)
866 if (!gmx_tng
->timePerFrameIsSet
&& gmx_tng
->lastTimeDataIsValid
&& gmx_tng
->lastStepDataIsValid
)
868 double deltaTime
= elapsedSeconds
- gmx_tng
->lastTime
;
869 std::int64_t deltaStep
= step
- gmx_tng
->lastStep
;
870 tng_time_per_frame_set(tng
, deltaTime
/ deltaStep
);
871 gmx_tng
->timePerFrameIsSet
= true;
874 tng_num_particles_get(tng
, &nParticles
);
875 if (nAtoms
!= (int)nParticles
)
877 tng_implicit_num_particles_set(tng
, nAtoms
);
880 if (bUseLossyCompression
)
882 compression
= TNG_TNG_COMPRESSION
;
886 compression
= TNG_GZIP_COMPRESSION
;
889 /* The writing is done using write_data, which writes float or double
890 * depending on the GROMACS compilation. */
893 GMX_ASSERT(box
, "Need a non-NULL box if positions are written");
895 if (write_data(tng
, step
, elapsedSeconds
,
896 reinterpret_cast<const real
*>(x
),
897 3, TNG_TRAJ_POSITIONS
, "POSITIONS",
898 TNG_PARTICLE_BLOCK_DATA
,
899 compression
) != TNG_SUCCESS
)
901 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
907 if (write_data(tng
, step
, elapsedSeconds
,
908 reinterpret_cast<const real
*>(v
),
909 3, TNG_TRAJ_VELOCITIES
, "VELOCITIES",
910 TNG_PARTICLE_BLOCK_DATA
,
911 compression
) != TNG_SUCCESS
)
913 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
919 /* TNG-MF1 compression only compresses positions and velocities. Use lossless
920 * compression for forces regardless of output mode */
921 if (write_data(tng
, step
, elapsedSeconds
,
922 reinterpret_cast<const real
*>(f
),
923 3, TNG_TRAJ_FORCES
, "FORCES",
924 TNG_PARTICLE_BLOCK_DATA
,
925 TNG_GZIP_COMPRESSION
) != TNG_SUCCESS
)
927 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
931 /* TNG-MF1 compression only compresses positions and velocities. Use lossless
932 * compression for lambdas and box shape regardless of output mode */
933 if (write_data(tng
, step
, elapsedSeconds
,
934 reinterpret_cast<const real
*>(box
),
935 9, TNG_TRAJ_BOX_SHAPE
, "BOX SHAPE",
936 TNG_NON_PARTICLE_BLOCK_DATA
,
937 TNG_GZIP_COMPRESSION
) != TNG_SUCCESS
)
939 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
942 if (write_data(tng
, step
, elapsedSeconds
,
943 reinterpret_cast<const real
*>(&lambda
),
944 1, TNG_GMX_LAMBDA
, "LAMBDAS",
945 TNG_NON_PARTICLE_BLOCK_DATA
,
946 TNG_GZIP_COMPRESSION
) != TNG_SUCCESS
)
948 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
951 // Update the data in the wrapper
953 gmx_tng
->lastStepDataIsValid
= true;
954 gmx_tng
->lastStep
= step
;
955 gmx_tng
->lastTimeDataIsValid
= true;
956 gmx_tng
->lastTime
= elapsedSeconds
;
958 GMX_UNUSED_VALUE(gmx_tng
);
959 GMX_UNUSED_VALUE(bUseLossyCompression
);
960 GMX_UNUSED_VALUE(step
);
961 GMX_UNUSED_VALUE(elapsedPicoSeconds
);
962 GMX_UNUSED_VALUE(lambda
);
963 GMX_UNUSED_VALUE(box
);
964 GMX_UNUSED_VALUE(nAtoms
);
971 void fflush_tng(gmx_tng_trajectory_t gmx_tng
)
978 tng_frame_set_premature_write(gmx_tng
->tng
, TNG_USE_HASH
);
980 GMX_UNUSED_VALUE(gmx_tng
);
984 float gmx_tng_get_time_of_final_frame(gmx_tng_trajectory_t gmx_tng
)
990 tng_trajectory_t tng
= gmx_tng
->tng
;
992 tng_num_frames_get(tng
, &nFrames
);
993 tng_util_time_of_frame_get(tng
, nFrames
- 1, &time
);
998 GMX_UNUSED_VALUE(gmx_tng
);
1003 void gmx_prepare_tng_writing(const char *filename
,
1005 gmx_tng_trajectory_t
*gmx_tng_input
,
1006 gmx_tng_trajectory_t
*gmx_tng_output
,
1008 const gmx_mtop_t
*mtop
,
1010 const char *indexGroupName
)
1013 tng_trajectory_t
*input
= (gmx_tng_input
&& *gmx_tng_input
) ? &(*gmx_tng_input
)->tng
: nullptr;
1014 /* FIXME after 5.0: Currently only standard block types are read */
1015 const int defaultNumIds
= 5;
1016 static gmx_int64_t fallbackIds
[defaultNumIds
] =
1018 TNG_TRAJ_BOX_SHAPE
, TNG_TRAJ_POSITIONS
,
1019 TNG_TRAJ_VELOCITIES
, TNG_TRAJ_FORCES
,
1022 static char fallbackNames
[defaultNumIds
][32] =
1024 "BOX SHAPE", "POSITIONS", "VELOCITIES",
1028 typedef tng_function_status (*set_writing_interval_func_pointer
)(tng_trajectory_t
,
1036 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_double_set
;
1038 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_set
;
1041 gmx_tng_open(filename
, mode
, gmx_tng_output
);
1042 tng_trajectory_t
*output
= &(*gmx_tng_output
)->tng
;
1044 /* Do we have an input file in TNG format? If so, then there's
1045 more data we can copy over, rather than having to improvise. */
1046 if (gmx_tng_input
&& *gmx_tng_input
)
1048 /* Set parameters (compression, time per frame, molecule
1049 * information, number of frames per frame set and writing
1050 * intervals of positions, box shape and lambdas) of the
1051 * output tng container based on their respective values int
1052 * the input tng container */
1053 double time
, compression_precision
;
1054 gmx_int64_t n_frames_per_frame_set
, interval
= -1;
1056 tng_compression_precision_get(*input
, &compression_precision
);
1057 tng_compression_precision_set(*output
, compression_precision
);
1058 // TODO make this configurable in a future version
1059 char compression_type
= TNG_TNG_COMPRESSION
;
1061 tng_molecule_system_copy(*input
, *output
);
1063 tng_time_per_frame_get(*input
, &time
);
1064 tng_time_per_frame_set(*output
, time
);
1065 // Since we have copied the value from the input TNG we should not change it again
1066 (*gmx_tng_output
)->timePerFrameIsSet
= true;
1068 tng_num_frames_per_frame_set_get(*input
, &n_frames_per_frame_set
);
1069 tng_num_frames_per_frame_set_set(*output
, n_frames_per_frame_set
);
1071 for (int i
= 0; i
< defaultNumIds
; i
++)
1073 if (tng_data_get_stride_length(*input
, fallbackIds
[i
], -1, &interval
)
1076 switch (fallbackIds
[i
])
1078 case TNG_TRAJ_POSITIONS
:
1079 case TNG_TRAJ_VELOCITIES
:
1080 set_writing_interval(*output
, interval
, 3, fallbackIds
[i
],
1081 fallbackNames
[i
], TNG_PARTICLE_BLOCK_DATA
,
1084 case TNG_TRAJ_FORCES
:
1085 set_writing_interval(*output
, interval
, 3, fallbackIds
[i
],
1086 fallbackNames
[i
], TNG_PARTICLE_BLOCK_DATA
,
1087 TNG_GZIP_COMPRESSION
);
1089 case TNG_TRAJ_BOX_SHAPE
:
1090 set_writing_interval(*output
, interval
, 9, fallbackIds
[i
],
1091 fallbackNames
[i
], TNG_NON_PARTICLE_BLOCK_DATA
,
1092 TNG_GZIP_COMPRESSION
);
1094 case TNG_GMX_LAMBDA
:
1095 set_writing_interval(*output
, interval
, 1, fallbackIds
[i
],
1096 fallbackNames
[i
], TNG_NON_PARTICLE_BLOCK_DATA
,
1097 TNG_GZIP_COMPRESSION
);
1108 /* TODO after trjconv is modularized: fix this so the user can
1109 change precision when they are doing an operation where
1110 this makes sense, and not otherwise.
1112 char compression = bUseLossyCompression ? TNG_TNG_COMPRESSION : TNG_GZIP_COMPRESSION;
1113 gmx_tng_set_compression_precision(*output, ndec2prec(nDecimalsOfPrecision));
1115 gmx_tng_add_mtop(*gmx_tng_output
, mtop
);
1116 tng_num_frames_per_frame_set_set(*output
, 1);
1119 if (index
&& nAtoms
> 0)
1121 gmx_tng_setup_atom_subgroup(*gmx_tng_output
, nAtoms
, index
, indexGroupName
);
1124 /* If for some reason there are more requested atoms than there are atoms in the
1125 * molecular system create a number of implicit atoms (without atom data) to
1126 * compensate for that. */
1129 tng_implicit_num_particles_set(*output
, nAtoms
);
1132 GMX_UNUSED_VALUE(filename
);
1133 GMX_UNUSED_VALUE(mode
);
1134 GMX_UNUSED_VALUE(gmx_tng_input
);
1135 GMX_UNUSED_VALUE(gmx_tng_output
);
1136 GMX_UNUSED_VALUE(nAtoms
);
1137 GMX_UNUSED_VALUE(mtop
);
1138 GMX_UNUSED_VALUE(index
);
1139 GMX_UNUSED_VALUE(indexGroupName
);
1143 void gmx_write_tng_from_trxframe(gmx_tng_trajectory_t gmx_tng_output
,
1144 const t_trxframe
*frame
,
1150 natoms
= frame
->natoms
;
1152 gmx_fwrite_tng(gmx_tng_output
,
1163 GMX_UNUSED_VALUE(gmx_tng_output
);
1164 GMX_UNUSED_VALUE(frame
);
1165 GMX_UNUSED_VALUE(natoms
);
1174 convert_array_to_real_array(void *from
,
1179 const char datatype
)
1183 const bool useDouble
= GMX_DOUBLE
;
1186 case TNG_FLOAT_DATA
:
1191 memcpy(to
, from
, nValues
* sizeof(real
) * nAtoms
);
1195 for (i
= 0; i
< nAtoms
; i
++)
1197 for (j
= 0; j
< nValues
; j
++)
1199 to
[i
*nValues
+j
] = reinterpret_cast<float *>(from
)[i
*nValues
+j
] * fact
;
1206 for (i
= 0; i
< nAtoms
; i
++)
1208 for (j
= 0; j
< nValues
; j
++)
1210 to
[i
*nValues
+j
] = reinterpret_cast<float *>(from
)[i
*nValues
+j
] * fact
;
1216 for (i
= 0; i
< nAtoms
; i
++)
1218 for (j
= 0; j
< nValues
; j
++)
1220 to
[i
*nValues
+j
] = reinterpret_cast<gmx_int64_t
*>(from
)[i
*nValues
+j
] * fact
;
1224 case TNG_DOUBLE_DATA
:
1225 if (sizeof(real
) == sizeof(double))
1229 memcpy(to
, from
, nValues
* sizeof(real
) * nAtoms
);
1233 for (i
= 0; i
< nAtoms
; i
++)
1235 for (j
= 0; j
< nValues
; j
++)
1237 to
[i
*nValues
+j
] = reinterpret_cast<double *>(from
)[i
*nValues
+j
] * fact
;
1244 for (i
= 0; i
< nAtoms
; i
++)
1246 for (j
= 0; j
< nValues
; j
++)
1248 to
[i
*nValues
+j
] = reinterpret_cast<double *>(from
)[i
*nValues
+j
] * fact
;
1254 gmx_incons("Illegal datatype when converting values to a real array!");
1260 real
getDistanceScaleFactor(gmx_tng_trajectory_t in
)
1262 gmx_int64_t exp
= -1;
1263 real distanceScaleFactor
;
1265 // TODO Hopefully, TNG 2.0 will do this kind of thing for us
1266 tng_distance_unit_exponential_get(in
->tng
, &exp
);
1268 // GROMACS expects distances in nm
1272 distanceScaleFactor
= NANO
/NANO
;
1275 distanceScaleFactor
= NANO
/ANGSTROM
;
1278 distanceScaleFactor
= pow(10.0, exp
+ 9.0);
1281 return distanceScaleFactor
;
1287 void gmx_tng_setup_atom_subgroup(gmx_tng_trajectory_t gmx_tng
,
1293 gmx_int64_t nAtoms
, cnt
, nMols
;
1294 tng_molecule_t mol
, iterMol
;
1298 tng_function_status stat
;
1299 tng_trajectory_t tng
= gmx_tng
->tng
;
1301 tng_num_particles_get(tng
, &nAtoms
);
1308 stat
= tng_molecule_find(tng
, name
, -1, &mol
);
1309 if (stat
== TNG_SUCCESS
)
1311 tng_molecule_num_atoms_get(tng
, mol
, &nAtoms
);
1312 tng_molecule_cnt_get(tng
, mol
, &cnt
);
1322 if (stat
== TNG_FAILURE
)
1324 /* The indexed atoms are added to one separate molecule. */
1325 tng_molecule_alloc(tng
, &mol
);
1326 tng_molecule_name_set(tng
, mol
, name
);
1327 tng_molecule_chain_add(tng
, mol
, "", &chain
);
1329 for (int i
= 0; i
< nind
; i
++)
1331 char temp_name
[256], temp_type
[256];
1333 /* Try to retrieve the residue name of the atom */
1334 stat
= tng_residue_name_of_particle_nr_get(tng
, ind
[i
], temp_name
, 256);
1335 if (stat
!= TNG_SUCCESS
)
1337 temp_name
[0] = '\0';
1339 /* Check if the molecule of the selection already contains this residue */
1340 if (tng_chain_residue_find(tng
, chain
, temp_name
, -1, &res
)
1343 tng_chain_residue_add(tng
, chain
, temp_name
, &res
);
1345 /* Try to find the original name and type of the atom */
1346 stat
= tng_atom_name_of_particle_nr_get(tng
, ind
[i
], temp_name
, 256);
1347 if (stat
!= TNG_SUCCESS
)
1349 temp_name
[0] = '\0';
1351 stat
= tng_atom_type_of_particle_nr_get(tng
, ind
[i
], temp_type
, 256);
1352 if (stat
!= TNG_SUCCESS
)
1354 temp_type
[0] = '\0';
1356 tng_residue_atom_w_id_add(tng
, res
, temp_name
, temp_type
, ind
[i
], &atom
);
1358 tng_molecule_existing_add(tng
, &mol
);
1360 /* Set the count of the molecule containing the selected atoms to 1 and all
1361 * other molecules to 0 */
1362 tng_molecule_cnt_set(tng
, mol
, 1);
1363 tng_num_molecule_types_get(tng
, &nMols
);
1364 for (gmx_int64_t k
= 0; k
< nMols
; k
++)
1366 tng_molecule_of_index_get(tng
, k
, &iterMol
);
1371 tng_molecule_cnt_set(tng
, iterMol
, 0);
1374 GMX_UNUSED_VALUE(gmx_tng
);
1375 GMX_UNUSED_VALUE(nind
);
1376 GMX_UNUSED_VALUE(ind
);
1377 GMX_UNUSED_VALUE(name
);
1381 /* TODO: If/when TNG acquires the ability to copy data blocks without
1382 * uncompressing them, then this implemenation should be reconsidered.
1383 * Ideally, gmx trjconv -f a.tng -o b.tng -b 10 -e 20 would be fast
1384 * and lose no information. */
1385 gmx_bool
gmx_read_next_tng_frame(gmx_tng_trajectory_t gmx_tng_input
,
1387 gmx_int64_t
*requestedIds
,
1388 int numRequestedIds
)
1391 tng_trajectory_t input
= gmx_tng_input
->tng
;
1392 gmx_bool bOK
= TRUE
;
1393 tng_function_status stat
;
1394 gmx_int64_t numberOfAtoms
= -1, frameNumber
= -1;
1395 gmx_int64_t nBlocks
, blockId
, *blockIds
= nullptr, codecId
;
1397 void *values
= nullptr;
1398 double frameTime
= -1.0;
1399 int size
, blockDependency
;
1401 const int defaultNumIds
= 5;
1402 static gmx_int64_t fallbackRequestedIds
[defaultNumIds
] =
1404 TNG_TRAJ_BOX_SHAPE
, TNG_TRAJ_POSITIONS
,
1405 TNG_TRAJ_VELOCITIES
, TNG_TRAJ_FORCES
,
1412 fr
->bLambda
= FALSE
;
1420 /* If no specific IDs were requested read all block types that can
1421 * currently be interpreted */
1422 if (!requestedIds
|| numRequestedIds
== 0)
1424 numRequestedIds
= defaultNumIds
;
1425 requestedIds
= fallbackRequestedIds
;
1428 stat
= tng_num_particles_get(input
, &numberOfAtoms
);
1429 if (stat
!= TNG_SUCCESS
)
1431 gmx_file("Cannot determine number of atoms from TNG file.");
1433 fr
->natoms
= numberOfAtoms
;
1435 bool nextFrameExists
= gmx_get_tng_data_block_types_of_next_frame(gmx_tng_input
,
1442 gmx::unique_cptr
<gmx_int64_t
, gmx::free_wrapper
> blockIdsGuard(blockIds
);
1443 if (!nextFrameExists
)
1453 for (gmx_int64_t i
= 0; i
< nBlocks
; i
++)
1455 blockId
= blockIds
[i
];
1456 tng_data_block_dependency_get(input
, blockId
, &blockDependency
);
1457 if (blockDependency
& TNG_PARTICLE_DEPENDENT
)
1459 stat
= tng_util_particle_data_next_frame_read(input
,
1468 stat
= tng_util_non_particle_data_next_frame_read(input
,
1475 if (stat
== TNG_CRITICAL
)
1477 gmx_file("Cannot read positions from TNG file.");
1480 else if (stat
== TNG_FAILURE
)
1486 case TNG_TRAJ_BOX_SHAPE
:
1490 size
= sizeof(gmx_int64_t
);
1492 case TNG_FLOAT_DATA
:
1493 size
= sizeof(float);
1495 case TNG_DOUBLE_DATA
:
1496 size
= sizeof(double);
1499 gmx_incons("Illegal datatype of box shape values!");
1501 for (int i
= 0; i
< DIM
; i
++)
1503 convert_array_to_real_array(reinterpret_cast<char *>(values
) + size
* i
* DIM
,
1504 reinterpret_cast<real
*>(fr
->box
[i
]),
1505 getDistanceScaleFactor(gmx_tng_input
),
1512 case TNG_TRAJ_POSITIONS
:
1513 srenew(fr
->x
, fr
->natoms
);
1514 convert_array_to_real_array(values
,
1515 reinterpret_cast<real
*>(fr
->x
),
1516 getDistanceScaleFactor(gmx_tng_input
),
1521 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &prec
);
1522 /* This must be updated if/when more lossy compression methods are added */
1523 if (codecId
== TNG_TNG_COMPRESSION
)
1529 case TNG_TRAJ_VELOCITIES
:
1530 srenew(fr
->v
, fr
->natoms
);
1531 convert_array_to_real_array(values
,
1533 getDistanceScaleFactor(gmx_tng_input
),
1538 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &prec
);
1539 /* This must be updated if/when more lossy compression methods are added */
1540 if (codecId
== TNG_TNG_COMPRESSION
)
1546 case TNG_TRAJ_FORCES
:
1547 srenew(fr
->f
, fr
->natoms
);
1548 convert_array_to_real_array(values
,
1549 reinterpret_cast<real
*>(fr
->f
),
1550 getDistanceScaleFactor(gmx_tng_input
),
1556 case TNG_GMX_LAMBDA
:
1559 case TNG_FLOAT_DATA
:
1560 fr
->lambda
= *(reinterpret_cast<float *>(values
));
1562 case TNG_DOUBLE_DATA
:
1563 fr
->lambda
= *(reinterpret_cast<double *>(values
));
1566 gmx_incons("Illegal datatype lambda value!");
1571 gmx_warning("Illegal block type! Currently GROMACS tools can only handle certain data types. Skipping block.");
1573 /* values does not have to be freed before reading next frame. It will
1574 * be reallocated if it is not NULL. */
1577 fr
->step
= frameNumber
;
1580 // Convert the time to ps
1581 fr
->time
= frameTime
/ PICO
;
1582 fr
->bTime
= (frameTime
> 0);
1584 // TODO This does not leak, but is not exception safe.
1585 /* values must be freed before leaving this function */
1590 GMX_UNUSED_VALUE(gmx_tng_input
);
1591 GMX_UNUSED_VALUE(fr
);
1592 GMX_UNUSED_VALUE(requestedIds
);
1593 GMX_UNUSED_VALUE(numRequestedIds
);
1598 void gmx_print_tng_molecule_system(gmx_tng_trajectory_t gmx_tng_input
,
1602 gmx_int64_t nMolecules
, nChains
, nResidues
, nAtoms
, nFramesRead
;
1603 gmx_int64_t strideLength
, nParticlesRead
, nValuesPerFrameRead
, *molCntList
;
1604 tng_molecule_t molecule
;
1606 tng_residue_t residue
;
1608 tng_function_status stat
;
1612 void *data
= nullptr;
1613 std::vector
<real
> atomCharges
;
1614 std::vector
<real
> atomMasses
;
1615 tng_trajectory_t input
= gmx_tng_input
->tng
;
1617 tng_num_molecule_types_get(input
, &nMolecules
);
1618 tng_molecule_cnt_list_get(input
, &molCntList
);
1619 /* Can the number of particles change in the trajectory or is it constant? */
1620 tng_num_particles_variable_get(input
, &varNAtoms
);
1622 for (gmx_int64_t i
= 0; i
< nMolecules
; i
++)
1624 tng_molecule_of_index_get(input
, i
, &molecule
);
1625 tng_molecule_name_get(input
, molecule
, str
, 256);
1626 if (varNAtoms
== TNG_CONSTANT_N_ATOMS
)
1628 if ((int)molCntList
[i
] == 0)
1632 fprintf(stream
, "Molecule: %s, count: %d\n", str
, (int)molCntList
[i
]);
1636 fprintf(stream
, "Molecule: %s\n", str
);
1638 tng_molecule_num_chains_get(input
, molecule
, &nChains
);
1641 for (gmx_int64_t j
= 0; j
< nChains
; j
++)
1643 tng_molecule_chain_of_index_get(input
, molecule
, j
, &chain
);
1644 tng_chain_name_get(input
, chain
, str
, 256);
1645 fprintf(stream
, "\tChain: %s\n", str
);
1646 tng_chain_num_residues_get(input
, chain
, &nResidues
);
1647 for (gmx_int64_t k
= 0; k
< nResidues
; k
++)
1649 tng_chain_residue_of_index_get(input
, chain
, k
, &residue
);
1650 tng_residue_name_get(input
, residue
, str
, 256);
1651 fprintf(stream
, "\t\tResidue: %s\n", str
);
1652 tng_residue_num_atoms_get(input
, residue
, &nAtoms
);
1653 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
1655 tng_residue_atom_of_index_get(input
, residue
, l
, &atom
);
1656 tng_atom_name_get(input
, atom
, str
, 256);
1657 fprintf(stream
, "\t\t\tAtom: %s", str
);
1658 tng_atom_type_get(input
, atom
, str
, 256);
1659 fprintf(stream
, " (%s)\n", str
);
1664 /* It is possible to have a molecule without chains, in which case
1665 * residues in the molecule can be iterated through without going
1666 * through chains. */
1669 tng_molecule_num_residues_get(input
, molecule
, &nResidues
);
1672 for (gmx_int64_t k
= 0; k
< nResidues
; k
++)
1674 tng_molecule_residue_of_index_get(input
, molecule
, k
, &residue
);
1675 tng_residue_name_get(input
, residue
, str
, 256);
1676 fprintf(stream
, "\t\tResidue: %s\n", str
);
1677 tng_residue_num_atoms_get(input
, residue
, &nAtoms
);
1678 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
1680 tng_residue_atom_of_index_get(input
, residue
, l
, &atom
);
1681 tng_atom_name_get(input
, atom
, str
, 256);
1682 fprintf(stream
, "\t\t\tAtom: %s", str
);
1683 tng_atom_type_get(input
, atom
, str
, 256);
1684 fprintf(stream
, " (%s)\n", str
);
1690 tng_molecule_num_atoms_get(input
, molecule
, &nAtoms
);
1691 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
1693 tng_molecule_atom_of_index_get(input
, molecule
, l
, &atom
);
1694 tng_atom_name_get(input
, atom
, str
, 256);
1695 fprintf(stream
, "\t\t\tAtom: %s", str
);
1696 tng_atom_type_get(input
, atom
, str
, 256);
1697 fprintf(stream
, " (%s)\n", str
);
1703 tng_num_particles_get(input
, &nAtoms
);
1704 stat
= tng_particle_data_vector_get(input
, TNG_TRAJ_PARTIAL_CHARGES
, &data
, &nFramesRead
,
1705 &strideLength
, &nParticlesRead
,
1706 &nValuesPerFrameRead
, &datatype
);
1707 if (stat
== TNG_SUCCESS
)
1709 atomCharges
.resize(nAtoms
);
1710 convert_array_to_real_array(data
,
1717 fprintf(stream
, "Atom Charges (%d):\n", int(nAtoms
));
1718 for (gmx_int64_t i
= 0; i
< nAtoms
; i
+= 10)
1720 fprintf(stream
, "Atom Charges [%8d-]=[", int(i
));
1721 for (gmx_int64_t j
= 0; (j
< 10 && i
+ j
< nAtoms
); j
++)
1723 fprintf(stream
, " %12.5e", atomCharges
[i
+ j
]);
1725 fprintf(stream
, "]\n");
1729 stat
= tng_particle_data_vector_get(input
, TNG_TRAJ_MASSES
, &data
, &nFramesRead
,
1730 &strideLength
, &nParticlesRead
,
1731 &nValuesPerFrameRead
, &datatype
);
1732 if (stat
== TNG_SUCCESS
)
1734 atomMasses
.resize(nAtoms
);
1735 convert_array_to_real_array(data
,
1742 fprintf(stream
, "Atom Masses (%d):\n", int(nAtoms
));
1743 for (gmx_int64_t i
= 0; i
< nAtoms
; i
+= 10)
1745 fprintf(stream
, "Atom Masses [%8d-]=[", int(i
));
1746 for (gmx_int64_t j
= 0; (j
< 10 && i
+ j
< nAtoms
); j
++)
1748 fprintf(stream
, " %12.5e", atomMasses
[i
+ j
]);
1750 fprintf(stream
, "]\n");
1756 GMX_UNUSED_VALUE(gmx_tng_input
);
1757 GMX_UNUSED_VALUE(stream
);
1761 gmx_bool
gmx_get_tng_data_block_types_of_next_frame(gmx_tng_trajectory_t gmx_tng_input
,
1764 gmx_int64_t
*requestedIds
,
1765 gmx_int64_t
*nextFrame
,
1766 gmx_int64_t
*nBlocks
,
1767 gmx_int64_t
**blockIds
)
1770 tng_function_status stat
;
1771 tng_trajectory_t input
= gmx_tng_input
->tng
;
1773 stat
= tng_util_trajectory_next_frame_present_data_blocks_find(input
, frame
,
1774 nRequestedIds
, requestedIds
,
1778 if (stat
== TNG_CRITICAL
)
1780 gmx_file("Cannot read TNG file. Cannot find data blocks of next frame.");
1782 else if (stat
== TNG_FAILURE
)
1788 GMX_UNUSED_VALUE(gmx_tng_input
);
1789 GMX_UNUSED_VALUE(frame
);
1790 GMX_UNUSED_VALUE(nRequestedIds
);
1791 GMX_UNUSED_VALUE(requestedIds
);
1792 GMX_UNUSED_VALUE(nextFrame
);
1793 GMX_UNUSED_VALUE(nBlocks
);
1794 GMX_UNUSED_VALUE(blockIds
);
1799 gmx_bool
gmx_get_tng_data_next_frame_of_block_type(gmx_tng_trajectory_t gmx_tng_input
,
1800 gmx_int64_t blockId
,
1802 gmx_int64_t
*frameNumber
,
1804 gmx_int64_t
*nValuesPerFrame
,
1805 gmx_int64_t
*nAtoms
,
1812 tng_function_status stat
;
1814 gmx_int64_t codecId
;
1815 int blockDependency
;
1816 void *data
= nullptr;
1818 tng_trajectory_t input
= gmx_tng_input
->tng
;
1820 stat
= tng_data_block_name_get(input
, blockId
, name
, maxLen
);
1821 if (stat
!= TNG_SUCCESS
)
1823 gmx_file("Cannot read next frame of TNG file");
1825 stat
= tng_data_block_dependency_get(input
, blockId
, &blockDependency
);
1826 if (stat
!= TNG_SUCCESS
)
1828 gmx_file("Cannot read next frame of TNG file");
1830 if (blockDependency
& TNG_PARTICLE_DEPENDENT
)
1832 tng_num_particles_get(input
, nAtoms
);
1833 stat
= tng_util_particle_data_next_frame_read(input
,
1842 *nAtoms
= 1; /* There are not actually any atoms, but it is used for
1843 allocating memory */
1844 stat
= tng_util_non_particle_data_next_frame_read(input
,
1851 if (stat
== TNG_CRITICAL
)
1853 gmx_file("Cannot read next frame of TNG file");
1855 if (stat
== TNG_FAILURE
)
1861 stat
= tng_data_block_num_values_per_frame_get(input
, blockId
, nValuesPerFrame
);
1862 if (stat
!= TNG_SUCCESS
)
1864 gmx_file("Cannot read next frame of TNG file");
1866 srenew(*values
, sizeof(real
) * *nValuesPerFrame
* *nAtoms
);
1867 convert_array_to_real_array(data
,
1869 getDistanceScaleFactor(gmx_tng_input
),
1874 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &localPrec
);
1876 /* This must be updated if/when more lossy compression methods are added */
1877 if (codecId
!= TNG_TNG_COMPRESSION
)
1891 GMX_UNUSED_VALUE(gmx_tng_input
);
1892 GMX_UNUSED_VALUE(blockId
);
1893 GMX_UNUSED_VALUE(values
);
1894 GMX_UNUSED_VALUE(frameNumber
);
1895 GMX_UNUSED_VALUE(frameTime
);
1896 GMX_UNUSED_VALUE(nValuesPerFrame
);
1897 GMX_UNUSED_VALUE(nAtoms
);
1898 GMX_UNUSED_VALUE(prec
);
1899 GMX_UNUSED_VALUE(name
);
1900 GMX_UNUSED_VALUE(maxLen
);
1901 GMX_UNUSED_VALUE(bOK
);