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"
69 using tng_trajectory_t
= void *;
72 /*! \brief Gromacs Wrapper around tng datatype
74 * This could in principle hold any GROMACS-specific requirements not yet
75 * implemented in or not relevant to the TNG library itself. However, for now
76 * we only use it to handle some shortcomings we have discovered, where the TNG
77 * API itself is a bit fragile and can end up overwriting data if called several
78 * times with the same frame number.
79 * The logic to determine the time per step was also a bit fragile. This is not
80 * critical, but since we anyway need a wrapper for ensuring unique frame
81 * numbers, we can also use it to store the time of the first step and use that
82 * to derive a slightly better/safer estimate of the time per step.
84 * At some future point where we have a second-generation TNG API we should
85 * consider removing this again.
87 struct gmx_tng_trajectory
89 tng_trajectory_t tng
; //!< Actual TNG handle (pointer)
90 bool lastStepDataIsValid
; //!< True if lastStep has been set
91 std::int64_t lastStep
; //!< Index/step used for last frame
92 bool lastTimeDataIsValid
; //!< True if lastTime has been set
93 double lastTime
; //!< Time of last frame (TNG unit is seconds)
94 bool timePerFrameIsSet
; //!< True if we have set the time per frame
98 static const char *modeToVerb(char mode
)
113 gmx_fatal(FARGS
, "Invalid file opening mode %c", mode
);
121 void gmx_tng_open(const char *filename
,
123 gmx_tng_trajectory_t
*gmx_tng
)
126 /* First check whether we have to make a backup,
127 * only for writing, not for read or append.
131 make_backup(filename
);
134 *gmx_tng
= new gmx_tng_trajectory
;
135 (*gmx_tng
)->lastStepDataIsValid
= false;
136 (*gmx_tng
)->lastTimeDataIsValid
= false;
137 (*gmx_tng
)->timePerFrameIsSet
= false;
138 tng_trajectory_t
* tng
= &(*gmx_tng
)->tng
;
140 /* tng must not be pointing at already allocated memory.
141 * Memory will be allocated by tng_util_trajectory_open() and must
142 * later on be freed by tng_util_trajectory_close(). */
143 if (TNG_SUCCESS
!= tng_util_trajectory_open(filename
, mode
, tng
))
145 /* TNG does return more than one degree of error, but there is
146 no use case for GROMACS handling the non-fatal errors
149 "File I/O error while opening %s for %s",
154 if (mode
== 'w' || mode
== 'a')
157 gmx_gethostname(hostname
, 256);
160 tng_first_computer_name_set(*tng
, hostname
);
164 tng_last_computer_name_set(*tng
, hostname
);
167 char programInfo
[256];
168 const char *precisionString
= "";
170 precisionString
= " (double precision)";
172 sprintf(programInfo
, "%.100s %.128s%.24s",
173 gmx::getProgramContext().displayName(),
174 gmx_version(), precisionString
);
177 tng_first_program_name_set(*tng
, programInfo
);
181 tng_last_program_name_set(*tng
, programInfo
);
185 if (!gmx_getusername(username
, 256))
189 tng_first_user_name_set(*tng
, username
);
193 tng_last_user_name_set(*tng
, username
);
194 tng_file_headers_write(*tng
, TNG_USE_HASH
);
199 gmx_file("GROMACS was compiled without TNG support, cannot handle this file type");
200 GMX_UNUSED_VALUE(filename
);
201 GMX_UNUSED_VALUE(mode
);
202 GMX_UNUSED_VALUE(gmx_tng
);
206 void gmx_tng_close(gmx_tng_trajectory_t
*gmx_tng
)
208 /* We have to check that tng is set because
209 * tng_util_trajectory_close wants to return a NULL in it, and
210 * gives a fatal error if it is NULL. */
212 if (gmx_tng
== nullptr || *gmx_tng
== nullptr)
216 tng_trajectory_t
* tng
= &(*gmx_tng
)->tng
;
220 tng_util_trajectory_close(tng
);
226 GMX_UNUSED_VALUE(gmx_tng
);
231 static void addTngMoleculeFromTopology(gmx_tng_trajectory_t gmx_tng
,
232 const char *moleculeName
,
233 const t_atoms
*atoms
,
234 gmx_int64_t numMolecules
,
235 tng_molecule_t
*tngMol
)
237 tng_trajectory_t tng
= gmx_tng
->tng
;
238 tng_chain_t tngChain
= nullptr;
239 tng_residue_t tngRes
= nullptr;
241 if (tng_molecule_add(tng
, moleculeName
, tngMol
) != TNG_SUCCESS
)
243 gmx_file("Cannot add molecule to TNG molecular system.");
246 for (int atomIndex
= 0; atomIndex
< atoms
->nr
; atomIndex
++)
248 const t_atom
*at
= &atoms
->atom
[atomIndex
];
249 /* FIXME: Currently the TNG API can only add atoms belonging to a
250 * residue and chain. Wait for TNG 2.0*/
253 const t_resinfo
*resInfo
= &atoms
->resinfo
[at
->resind
];
254 char chainName
[2] = {resInfo
->chainid
, 0};
255 tng_atom_t tngAtom
= nullptr;
260 prevAtom
= &atoms
->atom
[atomIndex
- 1];
267 /* If this is the first atom or if the residue changed add the
268 * residue to the TNG molecular system. */
269 if (!prevAtom
|| resInfo
!= &atoms
->resinfo
[prevAtom
->resind
])
271 /* If this is the first atom or if the chain changed add
272 * the chain to the TNG molecular system. */
273 if (!prevAtom
|| resInfo
->chainid
!=
274 atoms
->resinfo
[prevAtom
->resind
].chainid
)
276 tng_molecule_chain_add(tng
, *tngMol
, chainName
,
279 /* FIXME: When TNG supports both residue index and residue
280 * number the latter should be used. Wait for TNG 2.0*/
281 tng_chain_residue_add(tng
, tngChain
, *resInfo
->name
, &tngRes
);
283 tng_residue_atom_add(tng
, tngRes
, *(atoms
->atomname
[atomIndex
]), *(atoms
->atomtype
[atomIndex
]), &tngAtom
);
286 tng_molecule_cnt_set(tng
, *tngMol
, numMolecules
);
289 void gmx_tng_add_mtop(gmx_tng_trajectory_t gmx_tng
,
290 const gmx_mtop_t
*mtop
)
294 std::vector
<real
> atomCharges
;
295 std::vector
<real
> atomMasses
;
296 const t_ilist
*ilist
;
300 tng_trajectory_t tng
= gmx_tng
->tng
;
304 /* No topology information available to add. */
309 datatype
= TNG_DOUBLE_DATA
;
311 datatype
= TNG_FLOAT_DATA
;
314 atomCharges
.reserve(mtop
->natoms
);
315 atomMasses
.reserve(mtop
->natoms
);
317 for (int molIndex
= 0; molIndex
< mtop
->nmolblock
; molIndex
++)
319 tng_molecule_t tngMol
= nullptr;
320 const gmx_moltype_t
*molType
= &mtop
->moltype
[mtop
->molblock
[molIndex
].type
];
322 /* Add a molecule to the TNG trajectory with the same name as the
323 * current molecule. */
324 addTngMoleculeFromTopology(gmx_tng
,
327 mtop
->molblock
[molIndex
].nmol
,
330 /* Bonds have to be deduced from interactions (constraints etc). Different
331 * interactions have different sets of parameters. */
332 /* Constraints are specified using two atoms */
333 for (i
= 0; i
< F_NRE
; i
++)
337 ilist
= &molType
->ilist
[i
];
341 while (j
< ilist
->nr
)
343 tng_molecule_bond_add(tng
, tngMol
, ilist
->iatoms
[j
], ilist
->iatoms
[j
+1], &tngBond
);
349 /* Settle is described using three atoms */
350 ilist
= &molType
->ilist
[F_SETTLE
];
354 while (j
< ilist
->nr
)
356 tng_molecule_bond_add(tng
, tngMol
, ilist
->iatoms
[j
], ilist
->iatoms
[j
+1], &tngBond
);
357 tng_molecule_bond_add(tng
, tngMol
, ilist
->iatoms
[j
], ilist
->iatoms
[j
+2], &tngBond
);
361 /* First copy atom charges and masses, first atom by atom and then copy the memory for the molecule instances.
362 * FIXME: Atom B state data should also be written to TNG (v 2.0?) */
363 for (int atomCounter
= 0; atomCounter
< molType
->atoms
.nr
; atomCounter
++)
365 atomCharges
.push_back(molType
->atoms
.atom
[atomCounter
].q
);
366 atomMasses
.push_back(molType
->atoms
.atom
[atomCounter
].m
);
368 for (int molCounter
= 1; molCounter
< mtop
->molblock
[molIndex
].nmol
; molCounter
++)
370 std::copy_n(atomCharges
.end() - molType
->atoms
.nr
, molType
->atoms
.nr
, std::back_inserter(atomCharges
));
371 std::copy_n(atomMasses
.end() - molType
->atoms
.nr
, molType
->atoms
.nr
, std::back_inserter(atomMasses
));
374 /* Write the TNG data blocks. */
375 tng_particle_data_block_add(tng
, TNG_TRAJ_PARTIAL_CHARGES
, "PARTIAL CHARGES",
376 datatype
, TNG_NON_TRAJECTORY_BLOCK
,
377 1, 1, 1, 0, mtop
->natoms
,
378 TNG_GZIP_COMPRESSION
, atomCharges
.data());
379 tng_particle_data_block_add(tng
, TNG_TRAJ_MASSES
, "ATOM MASSES",
380 datatype
, TNG_NON_TRAJECTORY_BLOCK
,
381 1, 1, 1, 0, mtop
->natoms
,
382 TNG_GZIP_COMPRESSION
, atomMasses
.data());
385 /*! \libinternal \brief Compute greatest common divisor of n1 and n2
386 * if they are positive.
388 * If only one of n1 and n2 is positive, then return it.
389 * If neither n1 or n2 is positive, then return -1. */
391 greatest_common_divisor_if_positive(int n1
, int n2
)
395 return (0 >= n2
) ? -1 : n2
;
402 /* We have a non-trivial greatest common divisor to compute. */
403 return gmx_greatest_common_divisor(n1
, n2
);
406 /* By default try to write 100 frames (of actual output) in each frame set.
407 * This number is the number of outputs of the most frequently written data
408 * type per frame set.
409 * TODO for 5.1: Verify that 100 frames per frame set is efficient for most
410 * setups regarding compression efficiency and compression time. Make this
411 * a hidden command-line option? */
412 const int defaultFramesPerFrameSet
= 100;
414 /*! \libinternal \brief Set the number of frames per frame
415 * set according to output intervals.
416 * The default is that 100 frames are written of the data
417 * that is written most often. */
418 static void tng_set_frames_per_frame_set(gmx_tng_trajectory_t gmx_tng
,
419 const gmx_bool bUseLossyCompression
,
420 const t_inputrec
*ir
)
423 tng_trajectory_t tng
= gmx_tng
->tng
;
425 /* Set the number of frames per frame set to contain at least
426 * defaultFramesPerFrameSet of the lowest common denominator of
427 * the writing interval of positions and velocities. */
428 /* FIXME after 5.0: consider nstenergy also? */
429 if (bUseLossyCompression
)
431 gcd
= ir
->nstxout_compressed
;
435 gcd
= greatest_common_divisor_if_positive(ir
->nstxout
, ir
->nstvout
);
436 gcd
= greatest_common_divisor_if_positive(gcd
, ir
->nstfout
);
443 tng_num_frames_per_frame_set_set(tng
, gcd
* defaultFramesPerFrameSet
);
446 /*! \libinternal \brief Set the data-writing intervals, and number of
447 * frames per frame set */
448 static void set_writing_intervals(gmx_tng_trajectory_t gmx_tng
,
449 const gmx_bool bUseLossyCompression
,
450 const t_inputrec
*ir
)
452 tng_trajectory_t tng
= gmx_tng
->tng
;
454 /* Define pointers to specific writing functions depending on if we
455 * write float or double data */
456 typedef tng_function_status (*set_writing_interval_func_pointer
)(tng_trajectory_t
,
464 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_double_set
;
466 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_set
;
468 int xout
, vout
, fout
;
469 int gcd
= -1, lowest
= -1;
472 tng_set_frames_per_frame_set(gmx_tng
, bUseLossyCompression
, ir
);
474 if (bUseLossyCompression
)
476 xout
= ir
->nstxout_compressed
;
478 /* If there is no uncompressed coordinate output write forces
479 and velocities to the compressed tng file. */
490 compression
= TNG_TNG_COMPRESSION
;
497 compression
= TNG_GZIP_COMPRESSION
;
501 set_writing_interval(tng
, xout
, 3, TNG_TRAJ_POSITIONS
,
502 "POSITIONS", TNG_PARTICLE_BLOCK_DATA
,
504 /* TODO: if/when we write energies to TNG also, reconsider how
505 * and when box information is written, because GROMACS
506 * behaviour pre-5.0 was to write the box with every
507 * trajectory frame and every energy frame, and probably
508 * people depend on this. */
510 gcd
= greatest_common_divisor_if_positive(gcd
, xout
);
511 if (lowest
< 0 || xout
< lowest
)
518 set_writing_interval(tng
, vout
, 3, TNG_TRAJ_VELOCITIES
,
519 "VELOCITIES", TNG_PARTICLE_BLOCK_DATA
,
522 gcd
= greatest_common_divisor_if_positive(gcd
, vout
);
523 if (lowest
< 0 || vout
< lowest
)
530 set_writing_interval(tng
, fout
, 3, TNG_TRAJ_FORCES
,
531 "FORCES", TNG_PARTICLE_BLOCK_DATA
,
532 TNG_GZIP_COMPRESSION
);
534 gcd
= greatest_common_divisor_if_positive(gcd
, fout
);
535 if (lowest
< 0 || fout
< lowest
)
542 /* Lambdas and box shape written at an interval of the lowest common
543 denominator of other output */
544 set_writing_interval(tng
, gcd
, 1, TNG_GMX_LAMBDA
,
545 "LAMBDAS", TNG_NON_PARTICLE_BLOCK_DATA
,
546 TNG_GZIP_COMPRESSION
);
548 set_writing_interval(tng
, gcd
, 9, TNG_TRAJ_BOX_SHAPE
,
549 "BOX SHAPE", TNG_NON_PARTICLE_BLOCK_DATA
,
550 TNG_GZIP_COMPRESSION
);
551 if (gcd
< lowest
/ 10)
553 gmx_warning("The lowest common denominator of trajectory output is "
554 "every %d step(s), whereas the shortest output interval "
555 "is every %d steps.", gcd
, lowest
);
561 void gmx_tng_prepare_md_writing(gmx_tng_trajectory_t gmx_tng
,
562 const gmx_mtop_t
*mtop
,
563 const t_inputrec
*ir
)
566 gmx_tng_add_mtop(gmx_tng
, mtop
);
567 set_writing_intervals(gmx_tng
, FALSE
, ir
);
568 tng_time_per_frame_set(gmx_tng
->tng
, ir
->delta_t
* PICO
);
569 gmx_tng
->timePerFrameIsSet
= true;
571 GMX_UNUSED_VALUE(gmx_tng
);
572 GMX_UNUSED_VALUE(mtop
);
573 GMX_UNUSED_VALUE(ir
);
578 /* Check if all atoms in the molecule system are selected
579 * by a selection group of type specified by gtype. */
580 static gmx_bool
all_atoms_selected(const gmx_mtop_t
*mtop
,
583 const gmx_moltype_t
*molType
;
584 const t_atoms
*atoms
;
586 /* Iterate through all atoms in the molecule system and
587 * check if they belong to a selection group of the
589 for (int molIndex
= 0, i
= 0; molIndex
< mtop
->nmoltype
; molIndex
++)
591 molType
= &mtop
->moltype
[mtop
->molblock
[molIndex
].type
];
593 atoms
= &molType
->atoms
;
595 for (int j
= 0; j
< mtop
->molblock
[molIndex
].nmol
; j
++)
597 for (int atomIndex
= 0; atomIndex
< atoms
->nr
; atomIndex
++, i
++)
599 if (ggrpnr(&mtop
->groups
, gtype
, i
) != 0)
609 /* Create TNG molecules which will represent each of the selection
610 * groups for writing. But do that only if there is actually a
611 * specified selection group and it is not the whole system.
612 * TODO: Currently the only selection that is taken into account
613 * is egcCompressedX, but other selections should be added when
614 * e.g. writing energies is implemented.
616 static void add_selection_groups(gmx_tng_trajectory_t gmx_tng
,
617 const gmx_mtop_t
*mtop
)
619 const gmx_moltype_t
*molType
;
620 const t_atoms
*atoms
;
622 const t_resinfo
*resInfo
;
623 const t_ilist
*ilist
;
626 tng_molecule_t mol
, iterMol
;
633 tng_trajectory_t tng
= gmx_tng
->tng
;
635 /* TODO: When the TNG molecules block is more flexible TNG selection
636 * groups should not need all atoms specified. It should be possible
637 * just to specify what molecules are selected (and/or which atoms in
638 * the molecule) and how many (if applicable). */
640 /* If no atoms are selected we do not need to create a
641 * TNG selection group molecule. */
642 if (mtop
->groups
.ngrpnr
[egcCompressedX
] == 0)
647 /* If all atoms are selected we do not have to create a selection
648 * group molecule in the TNG molecule system. */
649 if (all_atoms_selected(mtop
, egcCompressedX
))
654 /* The name of the TNG molecule containing the selection group is the
655 * same as the name of the selection group. */
656 nameIndex
= *mtop
->groups
.grps
[egcCompressedX
].nm_ind
;
657 groupName
= *mtop
->groups
.grpname
[nameIndex
];
659 tng_molecule_alloc(tng
, &mol
);
660 tng_molecule_name_set(tng
, mol
, groupName
);
661 tng_molecule_chain_add(tng
, mol
, "", &chain
);
662 for (int molIndex
= 0, i
= 0; molIndex
< mtop
->nmolblock
; molIndex
++)
664 molType
= &mtop
->moltype
[mtop
->molblock
[molIndex
].type
];
666 atoms
= &molType
->atoms
;
668 for (int j
= 0; j
< mtop
->molblock
[molIndex
].nmol
; j
++)
670 bool bAtomsAdded
= FALSE
;
671 for (int atomIndex
= 0; atomIndex
< atoms
->nr
; atomIndex
++, i
++)
676 if (ggrpnr(&mtop
->groups
, egcCompressedX
, i
) != 0)
680 at
= &atoms
->atom
[atomIndex
];
683 resInfo
= &atoms
->resinfo
[at
->resind
];
684 /* FIXME: When TNG supports both residue index and residue
685 * number the latter should be used. */
686 res_name
= *resInfo
->name
;
687 res_id
= at
->resind
+ 1;
691 res_name
= (char *)"";
694 if (tng_chain_residue_find(tng
, chain
, res_name
, res_id
, &res
)
697 /* Since there is ONE chain for selection groups do not keep the
698 * original residue IDs - otherwise there might be conflicts. */
699 tng_chain_residue_add(tng
, chain
, res_name
, &res
);
701 tng_residue_atom_w_id_add(tng
, res
, *(atoms
->atomname
[atomIndex
]),
702 *(atoms
->atomtype
[atomIndex
]),
703 atom_offset
+ atomIndex
, &atom
);
709 for (int k
= 0; k
< F_NRE
; k
++)
713 ilist
= &molType
->ilist
[k
];
717 while (l
< ilist
->nr
)
720 atom1
= ilist
->iatoms
[l
] + atom_offset
;
721 atom2
= ilist
->iatoms
[l
+1] + atom_offset
;
722 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom1
) == 0 &&
723 ggrpnr(&mtop
->groups
, egcCompressedX
, atom2
) == 0)
725 tng_molecule_bond_add(tng
, mol
, ilist
->iatoms
[l
],
726 ilist
->iatoms
[l
+1], &tngBond
);
733 /* Settle is described using three atoms */
734 ilist
= &molType
->ilist
[F_SETTLE
];
738 while (l
< ilist
->nr
)
740 int atom1
, atom2
, atom3
;
741 atom1
= ilist
->iatoms
[l
] + atom_offset
;
742 atom2
= ilist
->iatoms
[l
+1] + atom_offset
;
743 atom3
= ilist
->iatoms
[l
+2] + atom_offset
;
744 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom1
) == 0)
746 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom2
) == 0)
748 tng_molecule_bond_add(tng
, mol
, atom1
,
751 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom3
) == 0)
753 tng_molecule_bond_add(tng
, mol
, atom1
,
761 atom_offset
+= atoms
->nr
;
764 tng_molecule_existing_add(tng
, &mol
);
765 tng_molecule_cnt_set(tng
, mol
, 1);
766 tng_num_molecule_types_get(tng
, &nMols
);
767 for (gmx_int64_t k
= 0; k
< nMols
; k
++)
769 tng_molecule_of_index_get(tng
, k
, &iterMol
);
774 tng_molecule_cnt_set(tng
, iterMol
, 0);
779 void gmx_tng_set_compression_precision(gmx_tng_trajectory_t gmx_tng
,
783 tng_compression_precision_set(gmx_tng
->tng
, prec
);
785 GMX_UNUSED_VALUE(gmx_tng
);
786 GMX_UNUSED_VALUE(prec
);
790 void gmx_tng_prepare_low_prec_writing(gmx_tng_trajectory_t gmx_tng
,
791 const gmx_mtop_t
*mtop
,
792 const t_inputrec
*ir
)
795 gmx_tng_add_mtop(gmx_tng
, mtop
);
796 add_selection_groups(gmx_tng
, mtop
);
797 set_writing_intervals(gmx_tng
, TRUE
, ir
);
798 tng_time_per_frame_set(gmx_tng
->tng
, ir
->delta_t
* PICO
);
799 gmx_tng
->timePerFrameIsSet
= true;
800 gmx_tng_set_compression_precision(gmx_tng
, ir
->x_compression_precision
);
802 GMX_UNUSED_VALUE(gmx_tng
);
803 GMX_UNUSED_VALUE(mtop
);
804 GMX_UNUSED_VALUE(ir
);
808 void gmx_fwrite_tng(gmx_tng_trajectory_t gmx_tng
,
809 const gmx_bool bUseLossyCompression
,
811 real elapsedPicoSeconds
,
820 typedef tng_function_status (*write_data_func_pointer
)(tng_trajectory_t
,
830 static write_data_func_pointer write_data
= tng_util_generic_with_time_double_write
;
832 static write_data_func_pointer write_data
= tng_util_generic_with_time_write
;
834 double elapsedSeconds
= elapsedPicoSeconds
* PICO
;
835 gmx_int64_t nParticles
;
841 /* This function might get called when the type of the
842 compressed trajectory is actually XTC. So we exit and move
846 tng_trajectory_t tng
= gmx_tng
->tng
;
848 // While the GROMACS interface to this routine specifies 'step', TNG itself
849 // only uses 'frame index' internally, although it suggests that it's a good
850 // idea to let this represent the step rather than just frame numbers.
852 // However, the way the GROMACS interface works, there are cases where
853 // the step is simply not set, so all frames rather have step=0.
854 // If we call the lower-level TNG routines with the same frame index
855 // (which is set from the step) multiple times, things will go very
856 // wrong (overwritten data).
858 // To avoid this, we require the step value to always be larger than the
859 // previous value, and if this is not true we simply set it to a value
860 // one unit larger than the previous step.
862 // Note: We do allow the user to specify any crazy value the want for the
863 // first step. If it is "not set", GROMACS will have used the value 0.
864 if (gmx_tng
->lastStepDataIsValid
&& step
<= gmx_tng
->lastStep
)
866 step
= gmx_tng
->lastStep
+ 1;
869 // Now that we have fixed the potentially incorrect step, we can also set
870 // the time per frame if it was not already done, and we have
871 // step/time information from the last frame so it is possible to calculate it.
872 // Note that we do allow the time to be the same for all steps, or even
873 // decreasing. In that case we will simply say the time per step is 0
874 // or negative. A bit strange, but a correct representation of the data :-)
875 if (!gmx_tng
->timePerFrameIsSet
&& gmx_tng
->lastTimeDataIsValid
&& gmx_tng
->lastStepDataIsValid
)
877 double deltaTime
= elapsedSeconds
- gmx_tng
->lastTime
;
878 std::int64_t deltaStep
= step
- gmx_tng
->lastStep
;
879 tng_time_per_frame_set(tng
, deltaTime
/ deltaStep
);
880 gmx_tng
->timePerFrameIsSet
= true;
883 tng_num_particles_get(tng
, &nParticles
);
884 if (nAtoms
!= (int)nParticles
)
886 tng_implicit_num_particles_set(tng
, nAtoms
);
889 if (bUseLossyCompression
)
891 compression
= TNG_TNG_COMPRESSION
;
895 compression
= TNG_GZIP_COMPRESSION
;
898 /* The writing is done using write_data, which writes float or double
899 * depending on the GROMACS compilation. */
902 GMX_ASSERT(box
, "Need a non-NULL box if positions are written");
904 if (write_data(tng
, step
, elapsedSeconds
,
905 reinterpret_cast<const real
*>(x
),
906 3, TNG_TRAJ_POSITIONS
, "POSITIONS",
907 TNG_PARTICLE_BLOCK_DATA
,
908 compression
) != TNG_SUCCESS
)
910 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
916 if (write_data(tng
, step
, elapsedSeconds
,
917 reinterpret_cast<const real
*>(v
),
918 3, TNG_TRAJ_VELOCITIES
, "VELOCITIES",
919 TNG_PARTICLE_BLOCK_DATA
,
920 compression
) != TNG_SUCCESS
)
922 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
928 /* TNG-MF1 compression only compresses positions and velocities. Use lossless
929 * compression for forces regardless of output mode */
930 if (write_data(tng
, step
, elapsedSeconds
,
931 reinterpret_cast<const real
*>(f
),
932 3, TNG_TRAJ_FORCES
, "FORCES",
933 TNG_PARTICLE_BLOCK_DATA
,
934 TNG_GZIP_COMPRESSION
) != TNG_SUCCESS
)
936 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
940 /* TNG-MF1 compression only compresses positions and velocities. Use lossless
941 * compression for lambdas and box shape regardless of output mode */
942 if (write_data(tng
, step
, elapsedSeconds
,
943 reinterpret_cast<const real
*>(box
),
944 9, TNG_TRAJ_BOX_SHAPE
, "BOX SHAPE",
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 if (write_data(tng
, step
, elapsedSeconds
,
952 reinterpret_cast<const real
*>(&lambda
),
953 1, TNG_GMX_LAMBDA
, "LAMBDAS",
954 TNG_NON_PARTICLE_BLOCK_DATA
,
955 TNG_GZIP_COMPRESSION
) != TNG_SUCCESS
)
957 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
960 // Update the data in the wrapper
962 gmx_tng
->lastStepDataIsValid
= true;
963 gmx_tng
->lastStep
= step
;
964 gmx_tng
->lastTimeDataIsValid
= true;
965 gmx_tng
->lastTime
= elapsedSeconds
;
967 GMX_UNUSED_VALUE(gmx_tng
);
968 GMX_UNUSED_VALUE(bUseLossyCompression
);
969 GMX_UNUSED_VALUE(step
);
970 GMX_UNUSED_VALUE(elapsedPicoSeconds
);
971 GMX_UNUSED_VALUE(lambda
);
972 GMX_UNUSED_VALUE(box
);
973 GMX_UNUSED_VALUE(nAtoms
);
980 void fflush_tng(gmx_tng_trajectory_t gmx_tng
)
987 tng_frame_set_premature_write(gmx_tng
->tng
, TNG_USE_HASH
);
989 GMX_UNUSED_VALUE(gmx_tng
);
993 float gmx_tng_get_time_of_final_frame(gmx_tng_trajectory_t gmx_tng
)
999 tng_trajectory_t tng
= gmx_tng
->tng
;
1001 tng_num_frames_get(tng
, &nFrames
);
1002 tng_util_time_of_frame_get(tng
, nFrames
- 1, &time
);
1004 fTime
= time
/ PICO
;
1007 GMX_UNUSED_VALUE(gmx_tng
);
1012 void gmx_prepare_tng_writing(const char *filename
,
1014 gmx_tng_trajectory_t
*gmx_tng_input
,
1015 gmx_tng_trajectory_t
*gmx_tng_output
,
1017 const gmx_mtop_t
*mtop
,
1019 const char *indexGroupName
)
1022 tng_trajectory_t
*input
= (gmx_tng_input
&& *gmx_tng_input
) ? &(*gmx_tng_input
)->tng
: nullptr;
1023 /* FIXME after 5.0: Currently only standard block types are read */
1024 const int defaultNumIds
= 5;
1025 static gmx_int64_t fallbackIds
[defaultNumIds
] =
1027 TNG_TRAJ_BOX_SHAPE
, TNG_TRAJ_POSITIONS
,
1028 TNG_TRAJ_VELOCITIES
, TNG_TRAJ_FORCES
,
1031 static char fallbackNames
[defaultNumIds
][32] =
1033 "BOX SHAPE", "POSITIONS", "VELOCITIES",
1037 typedef tng_function_status (*set_writing_interval_func_pointer
)(tng_trajectory_t
,
1045 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_double_set
;
1047 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_set
;
1050 gmx_tng_open(filename
, mode
, gmx_tng_output
);
1051 tng_trajectory_t
*output
= &(*gmx_tng_output
)->tng
;
1053 /* Do we have an input file in TNG format? If so, then there's
1054 more data we can copy over, rather than having to improvise. */
1055 if (gmx_tng_input
&& *gmx_tng_input
)
1057 /* Set parameters (compression, time per frame, molecule
1058 * information, number of frames per frame set and writing
1059 * intervals of positions, box shape and lambdas) of the
1060 * output tng container based on their respective values int
1061 * the input tng container */
1062 double time
, compression_precision
;
1063 gmx_int64_t n_frames_per_frame_set
, interval
= -1;
1065 tng_compression_precision_get(*input
, &compression_precision
);
1066 tng_compression_precision_set(*output
, compression_precision
);
1067 // TODO make this configurable in a future version
1068 char compression_type
= TNG_TNG_COMPRESSION
;
1070 tng_molecule_system_copy(*input
, *output
);
1072 tng_time_per_frame_get(*input
, &time
);
1073 tng_time_per_frame_set(*output
, time
);
1074 // Since we have copied the value from the input TNG we should not change it again
1075 (*gmx_tng_output
)->timePerFrameIsSet
= true;
1077 tng_num_frames_per_frame_set_get(*input
, &n_frames_per_frame_set
);
1078 tng_num_frames_per_frame_set_set(*output
, n_frames_per_frame_set
);
1080 for (int i
= 0; i
< defaultNumIds
; i
++)
1082 if (tng_data_get_stride_length(*input
, fallbackIds
[i
], -1, &interval
)
1085 switch (fallbackIds
[i
])
1087 case TNG_TRAJ_POSITIONS
:
1088 case TNG_TRAJ_VELOCITIES
:
1089 set_writing_interval(*output
, interval
, 3, fallbackIds
[i
],
1090 fallbackNames
[i
], TNG_PARTICLE_BLOCK_DATA
,
1093 case TNG_TRAJ_FORCES
:
1094 set_writing_interval(*output
, interval
, 3, fallbackIds
[i
],
1095 fallbackNames
[i
], TNG_PARTICLE_BLOCK_DATA
,
1096 TNG_GZIP_COMPRESSION
);
1098 case TNG_TRAJ_BOX_SHAPE
:
1099 set_writing_interval(*output
, interval
, 9, fallbackIds
[i
],
1100 fallbackNames
[i
], TNG_NON_PARTICLE_BLOCK_DATA
,
1101 TNG_GZIP_COMPRESSION
);
1103 case TNG_GMX_LAMBDA
:
1104 set_writing_interval(*output
, interval
, 1, fallbackIds
[i
],
1105 fallbackNames
[i
], TNG_NON_PARTICLE_BLOCK_DATA
,
1106 TNG_GZIP_COMPRESSION
);
1117 /* TODO after trjconv is modularized: fix this so the user can
1118 change precision when they are doing an operation where
1119 this makes sense, and not otherwise.
1121 char compression = bUseLossyCompression ? TNG_TNG_COMPRESSION : TNG_GZIP_COMPRESSION;
1122 gmx_tng_set_compression_precision(*output, ndec2prec(nDecimalsOfPrecision));
1124 gmx_tng_add_mtop(*gmx_tng_output
, mtop
);
1125 tng_num_frames_per_frame_set_set(*output
, 1);
1128 if (index
&& nAtoms
> 0)
1130 gmx_tng_setup_atom_subgroup(*gmx_tng_output
, nAtoms
, index
, indexGroupName
);
1133 /* If for some reason there are more requested atoms than there are atoms in the
1134 * molecular system create a number of implicit atoms (without atom data) to
1135 * compensate for that. */
1138 tng_implicit_num_particles_set(*output
, nAtoms
);
1141 GMX_UNUSED_VALUE(filename
);
1142 GMX_UNUSED_VALUE(mode
);
1143 GMX_UNUSED_VALUE(gmx_tng_input
);
1144 GMX_UNUSED_VALUE(gmx_tng_output
);
1145 GMX_UNUSED_VALUE(nAtoms
);
1146 GMX_UNUSED_VALUE(mtop
);
1147 GMX_UNUSED_VALUE(index
);
1148 GMX_UNUSED_VALUE(indexGroupName
);
1152 void gmx_write_tng_from_trxframe(gmx_tng_trajectory_t gmx_tng_output
,
1153 const t_trxframe
*frame
,
1159 natoms
= frame
->natoms
;
1161 gmx_fwrite_tng(gmx_tng_output
,
1172 GMX_UNUSED_VALUE(gmx_tng_output
);
1173 GMX_UNUSED_VALUE(frame
);
1174 GMX_UNUSED_VALUE(natoms
);
1183 convert_array_to_real_array(void *from
,
1188 const char datatype
)
1192 const bool useDouble
= GMX_DOUBLE
;
1195 case TNG_FLOAT_DATA
:
1200 memcpy(to
, from
, nValues
* sizeof(real
) * nAtoms
);
1204 for (i
= 0; i
< nAtoms
; i
++)
1206 for (j
= 0; j
< nValues
; j
++)
1208 to
[i
*nValues
+j
] = reinterpret_cast<float *>(from
)[i
*nValues
+j
] * fact
;
1215 for (i
= 0; i
< nAtoms
; i
++)
1217 for (j
= 0; j
< nValues
; j
++)
1219 to
[i
*nValues
+j
] = reinterpret_cast<float *>(from
)[i
*nValues
+j
] * fact
;
1225 for (i
= 0; i
< nAtoms
; i
++)
1227 for (j
= 0; j
< nValues
; j
++)
1229 to
[i
*nValues
+j
] = reinterpret_cast<gmx_int64_t
*>(from
)[i
*nValues
+j
] * fact
;
1233 case TNG_DOUBLE_DATA
:
1234 if (sizeof(real
) == sizeof(double))
1238 memcpy(to
, from
, nValues
* sizeof(real
) * nAtoms
);
1242 for (i
= 0; i
< nAtoms
; i
++)
1244 for (j
= 0; j
< nValues
; j
++)
1246 to
[i
*nValues
+j
] = reinterpret_cast<double *>(from
)[i
*nValues
+j
] * fact
;
1253 for (i
= 0; i
< nAtoms
; i
++)
1255 for (j
= 0; j
< nValues
; j
++)
1257 to
[i
*nValues
+j
] = reinterpret_cast<double *>(from
)[i
*nValues
+j
] * fact
;
1263 gmx_incons("Illegal datatype when converting values to a real array!");
1269 real
getDistanceScaleFactor(gmx_tng_trajectory_t in
)
1271 gmx_int64_t exp
= -1;
1272 real distanceScaleFactor
;
1274 // TODO Hopefully, TNG 2.0 will do this kind of thing for us
1275 tng_distance_unit_exponential_get(in
->tng
, &exp
);
1277 // GROMACS expects distances in nm
1281 distanceScaleFactor
= NANO
/NANO
;
1284 distanceScaleFactor
= NANO
/ANGSTROM
;
1287 distanceScaleFactor
= pow(10.0, exp
+ 9.0);
1290 return distanceScaleFactor
;
1296 void gmx_tng_setup_atom_subgroup(gmx_tng_trajectory_t gmx_tng
,
1302 gmx_int64_t nAtoms
, cnt
, nMols
;
1303 tng_molecule_t mol
, iterMol
;
1307 tng_function_status stat
;
1308 tng_trajectory_t tng
= gmx_tng
->tng
;
1310 tng_num_particles_get(tng
, &nAtoms
);
1317 stat
= tng_molecule_find(tng
, name
, -1, &mol
);
1318 if (stat
== TNG_SUCCESS
)
1320 tng_molecule_num_atoms_get(tng
, mol
, &nAtoms
);
1321 tng_molecule_cnt_get(tng
, mol
, &cnt
);
1331 if (stat
== TNG_FAILURE
)
1333 /* The indexed atoms are added to one separate molecule. */
1334 tng_molecule_alloc(tng
, &mol
);
1335 tng_molecule_name_set(tng
, mol
, name
);
1336 tng_molecule_chain_add(tng
, mol
, "", &chain
);
1338 for (int i
= 0; i
< nind
; i
++)
1340 char temp_name
[256], temp_type
[256];
1342 /* Try to retrieve the residue name of the atom */
1343 stat
= tng_residue_name_of_particle_nr_get(tng
, ind
[i
], temp_name
, 256);
1344 if (stat
!= TNG_SUCCESS
)
1346 temp_name
[0] = '\0';
1348 /* Check if the molecule of the selection already contains this residue */
1349 if (tng_chain_residue_find(tng
, chain
, temp_name
, -1, &res
)
1352 tng_chain_residue_add(tng
, chain
, temp_name
, &res
);
1354 /* Try to find the original name and type of the atom */
1355 stat
= tng_atom_name_of_particle_nr_get(tng
, ind
[i
], temp_name
, 256);
1356 if (stat
!= TNG_SUCCESS
)
1358 temp_name
[0] = '\0';
1360 stat
= tng_atom_type_of_particle_nr_get(tng
, ind
[i
], temp_type
, 256);
1361 if (stat
!= TNG_SUCCESS
)
1363 temp_type
[0] = '\0';
1365 tng_residue_atom_w_id_add(tng
, res
, temp_name
, temp_type
, ind
[i
], &atom
);
1367 tng_molecule_existing_add(tng
, &mol
);
1369 /* Set the count of the molecule containing the selected atoms to 1 and all
1370 * other molecules to 0 */
1371 tng_molecule_cnt_set(tng
, mol
, 1);
1372 tng_num_molecule_types_get(tng
, &nMols
);
1373 for (gmx_int64_t k
= 0; k
< nMols
; k
++)
1375 tng_molecule_of_index_get(tng
, k
, &iterMol
);
1380 tng_molecule_cnt_set(tng
, iterMol
, 0);
1383 GMX_UNUSED_VALUE(gmx_tng
);
1384 GMX_UNUSED_VALUE(nind
);
1385 GMX_UNUSED_VALUE(ind
);
1386 GMX_UNUSED_VALUE(name
);
1390 /* TODO: If/when TNG acquires the ability to copy data blocks without
1391 * uncompressing them, then this implemenation should be reconsidered.
1392 * Ideally, gmx trjconv -f a.tng -o b.tng -b 10 -e 20 would be fast
1393 * and lose no information. */
1394 gmx_bool
gmx_read_next_tng_frame(gmx_tng_trajectory_t gmx_tng_input
,
1396 gmx_int64_t
*requestedIds
,
1397 int numRequestedIds
)
1400 tng_trajectory_t input
= gmx_tng_input
->tng
;
1401 gmx_bool bOK
= TRUE
;
1402 tng_function_status stat
;
1403 gmx_int64_t numberOfAtoms
= -1, frameNumber
= -1;
1404 gmx_int64_t nBlocks
, blockId
, *blockIds
= nullptr, codecId
;
1406 void *values
= nullptr;
1407 double frameTime
= -1.0;
1408 int size
, blockDependency
;
1410 const int defaultNumIds
= 5;
1411 static gmx_int64_t fallbackRequestedIds
[defaultNumIds
] =
1413 TNG_TRAJ_BOX_SHAPE
, TNG_TRAJ_POSITIONS
,
1414 TNG_TRAJ_VELOCITIES
, TNG_TRAJ_FORCES
,
1421 fr
->bLambda
= FALSE
;
1429 /* If no specific IDs were requested read all block types that can
1430 * currently be interpreted */
1431 if (!requestedIds
|| numRequestedIds
== 0)
1433 numRequestedIds
= defaultNumIds
;
1434 requestedIds
= fallbackRequestedIds
;
1437 stat
= tng_num_particles_get(input
, &numberOfAtoms
);
1438 if (stat
!= TNG_SUCCESS
)
1440 gmx_file("Cannot determine number of atoms from TNG file.");
1442 fr
->natoms
= numberOfAtoms
;
1444 bool nextFrameExists
= gmx_get_tng_data_block_types_of_next_frame(gmx_tng_input
,
1451 gmx::unique_cptr
<gmx_int64_t
, gmx::free_wrapper
> blockIdsGuard(blockIds
);
1452 if (!nextFrameExists
)
1462 for (gmx_int64_t i
= 0; i
< nBlocks
; i
++)
1464 blockId
= blockIds
[i
];
1465 tng_data_block_dependency_get(input
, blockId
, &blockDependency
);
1466 if (blockDependency
& TNG_PARTICLE_DEPENDENT
)
1468 stat
= tng_util_particle_data_next_frame_read(input
,
1477 stat
= tng_util_non_particle_data_next_frame_read(input
,
1484 if (stat
== TNG_CRITICAL
)
1486 gmx_file("Cannot read positions from TNG file.");
1489 else if (stat
== TNG_FAILURE
)
1495 case TNG_TRAJ_BOX_SHAPE
:
1499 size
= sizeof(gmx_int64_t
);
1501 case TNG_FLOAT_DATA
:
1502 size
= sizeof(float);
1504 case TNG_DOUBLE_DATA
:
1505 size
= sizeof(double);
1508 gmx_incons("Illegal datatype of box shape values!");
1510 for (int i
= 0; i
< DIM
; i
++)
1512 convert_array_to_real_array(reinterpret_cast<char *>(values
) + size
* i
* DIM
,
1513 reinterpret_cast<real
*>(fr
->box
[i
]),
1514 getDistanceScaleFactor(gmx_tng_input
),
1521 case TNG_TRAJ_POSITIONS
:
1522 srenew(fr
->x
, fr
->natoms
);
1523 convert_array_to_real_array(values
,
1524 reinterpret_cast<real
*>(fr
->x
),
1525 getDistanceScaleFactor(gmx_tng_input
),
1530 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &prec
);
1531 /* This must be updated if/when more lossy compression methods are added */
1532 if (codecId
== TNG_TNG_COMPRESSION
)
1538 case TNG_TRAJ_VELOCITIES
:
1539 srenew(fr
->v
, fr
->natoms
);
1540 convert_array_to_real_array(values
,
1542 getDistanceScaleFactor(gmx_tng_input
),
1547 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &prec
);
1548 /* This must be updated if/when more lossy compression methods are added */
1549 if (codecId
== TNG_TNG_COMPRESSION
)
1555 case TNG_TRAJ_FORCES
:
1556 srenew(fr
->f
, fr
->natoms
);
1557 convert_array_to_real_array(values
,
1558 reinterpret_cast<real
*>(fr
->f
),
1559 getDistanceScaleFactor(gmx_tng_input
),
1565 case TNG_GMX_LAMBDA
:
1568 case TNG_FLOAT_DATA
:
1569 fr
->lambda
= *(reinterpret_cast<float *>(values
));
1571 case TNG_DOUBLE_DATA
:
1572 fr
->lambda
= *(reinterpret_cast<double *>(values
));
1575 gmx_incons("Illegal datatype lambda value!");
1580 gmx_warning("Illegal block type! Currently GROMACS tools can only handle certain data types. Skipping block.");
1582 /* values does not have to be freed before reading next frame. It will
1583 * be reallocated if it is not NULL. */
1586 fr
->step
= frameNumber
;
1589 // Convert the time to ps
1590 fr
->time
= frameTime
/ PICO
;
1591 fr
->bTime
= (frameTime
> 0);
1593 // TODO This does not leak, but is not exception safe.
1594 /* values must be freed before leaving this function */
1599 GMX_UNUSED_VALUE(gmx_tng_input
);
1600 GMX_UNUSED_VALUE(fr
);
1601 GMX_UNUSED_VALUE(requestedIds
);
1602 GMX_UNUSED_VALUE(numRequestedIds
);
1607 void gmx_print_tng_molecule_system(gmx_tng_trajectory_t gmx_tng_input
,
1611 gmx_int64_t nMolecules
, nChains
, nResidues
, nAtoms
, nFramesRead
;
1612 gmx_int64_t strideLength
, nParticlesRead
, nValuesPerFrameRead
, *molCntList
;
1613 tng_molecule_t molecule
;
1615 tng_residue_t residue
;
1617 tng_function_status stat
;
1621 void *data
= nullptr;
1622 std::vector
<real
> atomCharges
;
1623 std::vector
<real
> atomMasses
;
1624 tng_trajectory_t input
= gmx_tng_input
->tng
;
1626 tng_num_molecule_types_get(input
, &nMolecules
);
1627 tng_molecule_cnt_list_get(input
, &molCntList
);
1628 /* Can the number of particles change in the trajectory or is it constant? */
1629 tng_num_particles_variable_get(input
, &varNAtoms
);
1631 for (gmx_int64_t i
= 0; i
< nMolecules
; i
++)
1633 tng_molecule_of_index_get(input
, i
, &molecule
);
1634 tng_molecule_name_get(input
, molecule
, str
, 256);
1635 if (varNAtoms
== TNG_CONSTANT_N_ATOMS
)
1637 if ((int)molCntList
[i
] == 0)
1641 fprintf(stream
, "Molecule: %s, count: %d\n", str
, (int)molCntList
[i
]);
1645 fprintf(stream
, "Molecule: %s\n", str
);
1647 tng_molecule_num_chains_get(input
, molecule
, &nChains
);
1650 for (gmx_int64_t j
= 0; j
< nChains
; j
++)
1652 tng_molecule_chain_of_index_get(input
, molecule
, j
, &chain
);
1653 tng_chain_name_get(input
, chain
, str
, 256);
1654 fprintf(stream
, "\tChain: %s\n", str
);
1655 tng_chain_num_residues_get(input
, chain
, &nResidues
);
1656 for (gmx_int64_t k
= 0; k
< nResidues
; k
++)
1658 tng_chain_residue_of_index_get(input
, chain
, k
, &residue
);
1659 tng_residue_name_get(input
, residue
, str
, 256);
1660 fprintf(stream
, "\t\tResidue: %s\n", str
);
1661 tng_residue_num_atoms_get(input
, residue
, &nAtoms
);
1662 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
1664 tng_residue_atom_of_index_get(input
, residue
, l
, &atom
);
1665 tng_atom_name_get(input
, atom
, str
, 256);
1666 fprintf(stream
, "\t\t\tAtom: %s", str
);
1667 tng_atom_type_get(input
, atom
, str
, 256);
1668 fprintf(stream
, " (%s)\n", str
);
1673 /* It is possible to have a molecule without chains, in which case
1674 * residues in the molecule can be iterated through without going
1675 * through chains. */
1678 tng_molecule_num_residues_get(input
, molecule
, &nResidues
);
1681 for (gmx_int64_t k
= 0; k
< nResidues
; k
++)
1683 tng_molecule_residue_of_index_get(input
, molecule
, k
, &residue
);
1684 tng_residue_name_get(input
, residue
, str
, 256);
1685 fprintf(stream
, "\t\tResidue: %s\n", str
);
1686 tng_residue_num_atoms_get(input
, residue
, &nAtoms
);
1687 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
1689 tng_residue_atom_of_index_get(input
, residue
, l
, &atom
);
1690 tng_atom_name_get(input
, atom
, str
, 256);
1691 fprintf(stream
, "\t\t\tAtom: %s", str
);
1692 tng_atom_type_get(input
, atom
, str
, 256);
1693 fprintf(stream
, " (%s)\n", str
);
1699 tng_molecule_num_atoms_get(input
, molecule
, &nAtoms
);
1700 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
1702 tng_molecule_atom_of_index_get(input
, molecule
, l
, &atom
);
1703 tng_atom_name_get(input
, atom
, str
, 256);
1704 fprintf(stream
, "\t\t\tAtom: %s", str
);
1705 tng_atom_type_get(input
, atom
, str
, 256);
1706 fprintf(stream
, " (%s)\n", str
);
1712 tng_num_particles_get(input
, &nAtoms
);
1713 stat
= tng_particle_data_vector_get(input
, TNG_TRAJ_PARTIAL_CHARGES
, &data
, &nFramesRead
,
1714 &strideLength
, &nParticlesRead
,
1715 &nValuesPerFrameRead
, &datatype
);
1716 if (stat
== TNG_SUCCESS
)
1718 atomCharges
.resize(nAtoms
);
1719 convert_array_to_real_array(data
,
1726 fprintf(stream
, "Atom Charges (%d):\n", int(nAtoms
));
1727 for (gmx_int64_t i
= 0; i
< nAtoms
; i
+= 10)
1729 fprintf(stream
, "Atom Charges [%8d-]=[", int(i
));
1730 for (gmx_int64_t j
= 0; (j
< 10 && i
+ j
< nAtoms
); j
++)
1732 fprintf(stream
, " %12.5e", atomCharges
[i
+ j
]);
1734 fprintf(stream
, "]\n");
1738 stat
= tng_particle_data_vector_get(input
, TNG_TRAJ_MASSES
, &data
, &nFramesRead
,
1739 &strideLength
, &nParticlesRead
,
1740 &nValuesPerFrameRead
, &datatype
);
1741 if (stat
== TNG_SUCCESS
)
1743 atomMasses
.resize(nAtoms
);
1744 convert_array_to_real_array(data
,
1751 fprintf(stream
, "Atom Masses (%d):\n", int(nAtoms
));
1752 for (gmx_int64_t i
= 0; i
< nAtoms
; i
+= 10)
1754 fprintf(stream
, "Atom Masses [%8d-]=[", int(i
));
1755 for (gmx_int64_t j
= 0; (j
< 10 && i
+ j
< nAtoms
); j
++)
1757 fprintf(stream
, " %12.5e", atomMasses
[i
+ j
]);
1759 fprintf(stream
, "]\n");
1765 GMX_UNUSED_VALUE(gmx_tng_input
);
1766 GMX_UNUSED_VALUE(stream
);
1770 gmx_bool
gmx_get_tng_data_block_types_of_next_frame(gmx_tng_trajectory_t gmx_tng_input
,
1773 gmx_int64_t
*requestedIds
,
1774 gmx_int64_t
*nextFrame
,
1775 gmx_int64_t
*nBlocks
,
1776 gmx_int64_t
**blockIds
)
1779 tng_function_status stat
;
1780 tng_trajectory_t input
= gmx_tng_input
->tng
;
1782 stat
= tng_util_trajectory_next_frame_present_data_blocks_find(input
, frame
,
1783 nRequestedIds
, requestedIds
,
1787 if (stat
== TNG_CRITICAL
)
1789 gmx_file("Cannot read TNG file. Cannot find data blocks of next frame.");
1791 else if (stat
== TNG_FAILURE
)
1797 GMX_UNUSED_VALUE(gmx_tng_input
);
1798 GMX_UNUSED_VALUE(frame
);
1799 GMX_UNUSED_VALUE(nRequestedIds
);
1800 GMX_UNUSED_VALUE(requestedIds
);
1801 GMX_UNUSED_VALUE(nextFrame
);
1802 GMX_UNUSED_VALUE(nBlocks
);
1803 GMX_UNUSED_VALUE(blockIds
);
1808 gmx_bool
gmx_get_tng_data_next_frame_of_block_type(gmx_tng_trajectory_t gmx_tng_input
,
1809 gmx_int64_t blockId
,
1811 gmx_int64_t
*frameNumber
,
1813 gmx_int64_t
*nValuesPerFrame
,
1814 gmx_int64_t
*nAtoms
,
1821 tng_function_status stat
;
1823 gmx_int64_t codecId
;
1824 int blockDependency
;
1825 void *data
= nullptr;
1827 tng_trajectory_t input
= gmx_tng_input
->tng
;
1829 stat
= tng_data_block_name_get(input
, blockId
, name
, maxLen
);
1830 if (stat
!= TNG_SUCCESS
)
1832 gmx_file("Cannot read next frame of TNG file");
1834 stat
= tng_data_block_dependency_get(input
, blockId
, &blockDependency
);
1835 if (stat
!= TNG_SUCCESS
)
1837 gmx_file("Cannot read next frame of TNG file");
1839 if (blockDependency
& TNG_PARTICLE_DEPENDENT
)
1841 tng_num_particles_get(input
, nAtoms
);
1842 stat
= tng_util_particle_data_next_frame_read(input
,
1851 *nAtoms
= 1; /* There are not actually any atoms, but it is used for
1852 allocating memory */
1853 stat
= tng_util_non_particle_data_next_frame_read(input
,
1860 if (stat
== TNG_CRITICAL
)
1862 gmx_file("Cannot read next frame of TNG file");
1864 if (stat
== TNG_FAILURE
)
1870 stat
= tng_data_block_num_values_per_frame_get(input
, blockId
, nValuesPerFrame
);
1871 if (stat
!= TNG_SUCCESS
)
1873 gmx_file("Cannot read next frame of TNG file");
1875 srenew(*values
, sizeof(real
) * *nValuesPerFrame
* *nAtoms
);
1876 convert_array_to_real_array(data
,
1878 getDistanceScaleFactor(gmx_tng_input
),
1883 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &localPrec
);
1885 /* This must be updated if/when more lossy compression methods are added */
1886 if (codecId
!= TNG_TNG_COMPRESSION
)
1900 GMX_UNUSED_VALUE(gmx_tng_input
);
1901 GMX_UNUSED_VALUE(blockId
);
1902 GMX_UNUSED_VALUE(values
);
1903 GMX_UNUSED_VALUE(frameNumber
);
1904 GMX_UNUSED_VALUE(frameTime
);
1905 GMX_UNUSED_VALUE(nValuesPerFrame
);
1906 GMX_UNUSED_VALUE(nAtoms
);
1907 GMX_UNUSED_VALUE(prec
);
1908 GMX_UNUSED_VALUE(name
);
1909 GMX_UNUSED_VALUE(maxLen
);
1910 GMX_UNUSED_VALUE(bOK
);