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
91 int boxOutputInterval
; //!< Number of steps between the output of box size
92 int lambdaOutputInterval
; //!< Number of steps between the output of lambdas
95 static const char *modeToVerb(char mode
)
110 gmx_fatal(FARGS
, "Invalid file opening mode %c", mode
);
117 void gmx_tng_open(const char *filename
,
119 gmx_tng_trajectory_t
*gmx_tng
)
122 /* First check whether we have to make a backup,
123 * only for writing, not for read or append.
127 make_backup(filename
);
130 *gmx_tng
= new gmx_tng_trajectory
;
131 (*gmx_tng
)->lastStepDataIsValid
= false;
132 (*gmx_tng
)->lastTimeDataIsValid
= false;
133 (*gmx_tng
)->timePerFrameIsSet
= false;
134 tng_trajectory_t
* tng
= &(*gmx_tng
)->tng
;
136 /* tng must not be pointing at already allocated memory.
137 * Memory will be allocated by tng_util_trajectory_open() and must
138 * later on be freed by tng_util_trajectory_close(). */
139 if (TNG_SUCCESS
!= tng_util_trajectory_open(filename
, mode
, tng
))
141 /* TNG does return more than one degree of error, but there is
142 no use case for GROMACS handling the non-fatal errors
145 "File I/O error while opening %s for %s",
150 if (mode
== 'w' || mode
== 'a')
153 gmx_gethostname(hostname
, 256);
156 tng_first_computer_name_set(*tng
, hostname
);
160 tng_last_computer_name_set(*tng
, hostname
);
163 char programInfo
[256];
164 const char *precisionString
= "";
166 precisionString
= " (double precision)";
168 sprintf(programInfo
, "%.100s %.128s%.24s",
169 gmx::getProgramContext().displayName(),
170 gmx_version(), precisionString
);
173 tng_first_program_name_set(*tng
, programInfo
);
177 tng_last_program_name_set(*tng
, programInfo
);
181 if (!gmx_getusername(username
, 256))
185 tng_first_user_name_set(*tng
, username
);
189 tng_last_user_name_set(*tng
, username
);
190 tng_file_headers_write(*tng
, TNG_USE_HASH
);
195 gmx_file("GROMACS was compiled without TNG support, cannot handle this file type");
196 GMX_UNUSED_VALUE(filename
);
197 GMX_UNUSED_VALUE(mode
);
198 GMX_UNUSED_VALUE(gmx_tng
);
202 void gmx_tng_close(gmx_tng_trajectory_t
*gmx_tng
)
204 /* We have to check that tng is set because
205 * tng_util_trajectory_close wants to return a NULL in it, and
206 * gives a fatal error if it is NULL. */
208 if (gmx_tng
== nullptr || *gmx_tng
== nullptr)
212 tng_trajectory_t
* tng
= &(*gmx_tng
)->tng
;
216 tng_util_trajectory_close(tng
);
222 GMX_UNUSED_VALUE(gmx_tng
);
227 static void addTngMoleculeFromTopology(gmx_tng_trajectory_t gmx_tng
,
228 const char *moleculeName
,
229 const t_atoms
*atoms
,
230 gmx_int64_t numMolecules
,
231 tng_molecule_t
*tngMol
)
233 tng_trajectory_t tng
= gmx_tng
->tng
;
234 tng_chain_t tngChain
= nullptr;
235 tng_residue_t tngRes
= nullptr;
237 if (tng_molecule_add(tng
, moleculeName
, tngMol
) != TNG_SUCCESS
)
239 gmx_file("Cannot add molecule to TNG molecular system.");
242 for (int atomIndex
= 0; atomIndex
< atoms
->nr
; atomIndex
++)
244 const t_atom
*at
= &atoms
->atom
[atomIndex
];
245 /* FIXME: Currently the TNG API can only add atoms belonging to a
246 * residue and chain. Wait for TNG 2.0*/
249 const t_resinfo
*resInfo
= &atoms
->resinfo
[at
->resind
];
250 char chainName
[2] = {resInfo
->chainid
, 0};
251 tng_atom_t tngAtom
= nullptr;
256 prevAtom
= &atoms
->atom
[atomIndex
- 1];
263 /* If this is the first atom or if the residue changed add the
264 * residue to the TNG molecular system. */
265 if (!prevAtom
|| resInfo
!= &atoms
->resinfo
[prevAtom
->resind
])
267 /* If this is the first atom or if the chain changed add
268 * the chain to the TNG molecular system. */
269 if (!prevAtom
|| resInfo
->chainid
!=
270 atoms
->resinfo
[prevAtom
->resind
].chainid
)
272 tng_molecule_chain_add(tng
, *tngMol
, chainName
,
275 /* FIXME: When TNG supports both residue index and residue
276 * number the latter should be used. Wait for TNG 2.0*/
277 tng_chain_residue_add(tng
, tngChain
, *resInfo
->name
, &tngRes
);
279 tng_residue_atom_add(tng
, tngRes
, *(atoms
->atomname
[atomIndex
]), *(atoms
->atomtype
[atomIndex
]), &tngAtom
);
282 tng_molecule_cnt_set(tng
, *tngMol
, numMolecules
);
285 void gmx_tng_add_mtop(gmx_tng_trajectory_t gmx_tng
,
286 const gmx_mtop_t
*mtop
)
290 std::vector
<real
> atomCharges
;
291 std::vector
<real
> atomMasses
;
292 const t_ilist
*ilist
;
296 tng_trajectory_t tng
= gmx_tng
->tng
;
300 /* No topology information available to add. */
305 datatype
= TNG_DOUBLE_DATA
;
307 datatype
= TNG_FLOAT_DATA
;
310 atomCharges
.reserve(mtop
->natoms
);
311 atomMasses
.reserve(mtop
->natoms
);
313 for (const gmx_molblock_t
&molBlock
: mtop
->molblock
)
315 tng_molecule_t tngMol
= nullptr;
316 const gmx_moltype_t
*molType
= &mtop
->moltype
[molBlock
.type
];
318 /* Add a molecule to the TNG trajectory with the same name as the
319 * current molecule. */
320 addTngMoleculeFromTopology(gmx_tng
,
326 /* Bonds have to be deduced from interactions (constraints etc). Different
327 * interactions have different sets of parameters. */
328 /* Constraints are specified using two atoms */
329 for (i
= 0; i
< F_NRE
; i
++)
333 ilist
= &molType
->ilist
[i
];
337 while (j
< ilist
->nr
)
339 tng_molecule_bond_add(tng
, tngMol
, ilist
->iatoms
[j
], ilist
->iatoms
[j
+1], &tngBond
);
345 /* Settle is described using three atoms */
346 ilist
= &molType
->ilist
[F_SETTLE
];
350 while (j
< ilist
->nr
)
352 tng_molecule_bond_add(tng
, tngMol
, ilist
->iatoms
[j
], ilist
->iatoms
[j
+1], &tngBond
);
353 tng_molecule_bond_add(tng
, tngMol
, ilist
->iatoms
[j
], ilist
->iatoms
[j
+2], &tngBond
);
357 /* First copy atom charges and masses, first atom by atom and then copy the memory for the molecule instances.
358 * FIXME: Atom B state data should also be written to TNG (v 2.0?) */
359 for (int atomCounter
= 0; atomCounter
< molType
->atoms
.nr
; atomCounter
++)
361 atomCharges
.push_back(molType
->atoms
.atom
[atomCounter
].q
);
362 atomMasses
.push_back(molType
->atoms
.atom
[atomCounter
].m
);
364 for (int molCounter
= 1; molCounter
< molBlock
.nmol
; molCounter
++)
366 std::copy_n(atomCharges
.end() - molType
->atoms
.nr
, molType
->atoms
.nr
, std::back_inserter(atomCharges
));
367 std::copy_n(atomMasses
.end() - molType
->atoms
.nr
, molType
->atoms
.nr
, std::back_inserter(atomMasses
));
370 /* Write the TNG data blocks. */
371 tng_particle_data_block_add(tng
, TNG_TRAJ_PARTIAL_CHARGES
, "PARTIAL CHARGES",
372 datatype
, TNG_NON_TRAJECTORY_BLOCK
,
373 1, 1, 1, 0, mtop
->natoms
,
374 TNG_GZIP_COMPRESSION
, atomCharges
.data());
375 tng_particle_data_block_add(tng
, TNG_TRAJ_MASSES
, "ATOM MASSES",
376 datatype
, TNG_NON_TRAJECTORY_BLOCK
,
377 1, 1, 1, 0, mtop
->natoms
,
378 TNG_GZIP_COMPRESSION
, atomMasses
.data());
381 /*! \libinternal \brief Compute greatest common divisor of n1 and n2
382 * if they are positive.
384 * If only one of n1 and n2 is positive, then return it.
385 * If neither n1 or n2 is positive, then return -1. */
387 greatest_common_divisor_if_positive(int n1
, int n2
)
391 return (0 >= n2
) ? -1 : n2
;
398 /* We have a non-trivial greatest common divisor to compute. */
399 return gmx_greatest_common_divisor(n1
, n2
);
402 /* By default try to write 100 frames (of actual output) in each frame set.
403 * This number is the number of outputs of the most frequently written data
404 * type per frame set.
405 * TODO for 5.1: Verify that 100 frames per frame set is efficient for most
406 * setups regarding compression efficiency and compression time. Make this
407 * a hidden command-line option? */
408 const int defaultFramesPerFrameSet
= 100;
410 /*! \libinternal \brief Set the number of frames per frame
411 * set according to output intervals.
412 * The default is that 100 frames are written of the data
413 * that is written most often. */
414 static void tng_set_frames_per_frame_set(gmx_tng_trajectory_t gmx_tng
,
415 const gmx_bool bUseLossyCompression
,
416 const t_inputrec
*ir
)
419 tng_trajectory_t tng
= gmx_tng
->tng
;
421 /* Set the number of frames per frame set to contain at least
422 * defaultFramesPerFrameSet of the lowest common denominator of
423 * the writing interval of positions and velocities. */
424 /* FIXME after 5.0: consider nstenergy also? */
425 if (bUseLossyCompression
)
427 gcd
= ir
->nstxout_compressed
;
431 gcd
= greatest_common_divisor_if_positive(ir
->nstxout
, ir
->nstvout
);
432 gcd
= greatest_common_divisor_if_positive(gcd
, ir
->nstfout
);
439 tng_num_frames_per_frame_set_set(tng
, gcd
* defaultFramesPerFrameSet
);
442 /*! \libinternal \brief Set the data-writing intervals, and number of
443 * frames per frame set */
444 static void set_writing_intervals(gmx_tng_trajectory_t gmx_tng
,
445 const gmx_bool bUseLossyCompression
,
446 const t_inputrec
*ir
)
448 tng_trajectory_t tng
= gmx_tng
->tng
;
450 /* Define pointers to specific writing functions depending on if we
451 * write float or double data */
452 typedef tng_function_status (*set_writing_interval_func_pointer
)(tng_trajectory_t
,
460 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_double_set
;
462 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_set
;
464 int xout
, vout
, fout
;
465 int gcd
= -1, lowest
= -1;
468 tng_set_frames_per_frame_set(gmx_tng
, bUseLossyCompression
, ir
);
470 if (bUseLossyCompression
)
472 xout
= ir
->nstxout_compressed
;
474 /* If there is no uncompressed coordinate output write forces
475 and velocities to the compressed tng file. */
486 compression
= TNG_TNG_COMPRESSION
;
493 compression
= TNG_GZIP_COMPRESSION
;
497 set_writing_interval(tng
, xout
, 3, TNG_TRAJ_POSITIONS
,
498 "POSITIONS", TNG_PARTICLE_BLOCK_DATA
,
500 /* TODO: if/when we write energies to TNG also, reconsider how
501 * and when box information is written, because GROMACS
502 * behaviour pre-5.0 was to write the box with every
503 * trajectory frame and every energy frame, and probably
504 * people depend on this. */
506 gcd
= greatest_common_divisor_if_positive(gcd
, xout
);
507 if (lowest
< 0 || xout
< lowest
)
514 set_writing_interval(tng
, vout
, 3, TNG_TRAJ_VELOCITIES
,
515 "VELOCITIES", TNG_PARTICLE_BLOCK_DATA
,
518 gcd
= greatest_common_divisor_if_positive(gcd
, vout
);
519 if (lowest
< 0 || vout
< lowest
)
526 set_writing_interval(tng
, fout
, 3, TNG_TRAJ_FORCES
,
527 "FORCES", TNG_PARTICLE_BLOCK_DATA
,
528 TNG_GZIP_COMPRESSION
);
530 gcd
= greatest_common_divisor_if_positive(gcd
, fout
);
531 if (lowest
< 0 || fout
< lowest
)
538 /* Lambdas and box shape written at an interval of the lowest common
539 denominator of other output */
540 set_writing_interval(tng
, gcd
, 1, TNG_GMX_LAMBDA
,
541 "LAMBDAS", TNG_NON_PARTICLE_BLOCK_DATA
,
542 TNG_GZIP_COMPRESSION
);
544 set_writing_interval(tng
, gcd
, 9, TNG_TRAJ_BOX_SHAPE
,
545 "BOX SHAPE", TNG_NON_PARTICLE_BLOCK_DATA
,
546 TNG_GZIP_COMPRESSION
);
547 gmx_tng
->lambdaOutputInterval
= gcd
;
548 gmx_tng
->boxOutputInterval
= gcd
;
549 if (gcd
< lowest
/ 10)
551 gmx_warning("The lowest common denominator of trajectory output is "
552 "every %d step(s), whereas the shortest output interval "
553 "is every %d steps.", gcd
, lowest
);
559 void gmx_tng_prepare_md_writing(gmx_tng_trajectory_t gmx_tng
,
560 const gmx_mtop_t
*mtop
,
561 const t_inputrec
*ir
)
564 gmx_tng_add_mtop(gmx_tng
, mtop
);
565 set_writing_intervals(gmx_tng
, FALSE
, ir
);
566 tng_time_per_frame_set(gmx_tng
->tng
, ir
->delta_t
* PICO
);
567 gmx_tng
->timePerFrameIsSet
= true;
569 GMX_UNUSED_VALUE(gmx_tng
);
570 GMX_UNUSED_VALUE(mtop
);
571 GMX_UNUSED_VALUE(ir
);
576 /* Check if all atoms in the molecule system are selected
577 * by a selection group of type specified by gtype. */
578 static gmx_bool
all_atoms_selected(const gmx_mtop_t
*mtop
,
581 /* Iterate through all atoms in the molecule system and
582 * check if they belong to a selection group of the
585 for (const gmx_molblock_t
&molBlock
: mtop
->molblock
)
587 const gmx_moltype_t
&molType
= mtop
->moltype
[molBlock
.type
];
588 const t_atoms
&atoms
= molType
.atoms
;
590 for (int j
= 0; j
< molBlock
.nmol
; j
++)
592 for (int atomIndex
= 0; atomIndex
< atoms
.nr
; atomIndex
++, i
++)
594 if (ggrpnr(&mtop
->groups
, gtype
, i
) != 0)
604 /* Create TNG molecules which will represent each of the selection
605 * groups for writing. But do that only if there is actually a
606 * specified selection group and it is not the whole system.
607 * TODO: Currently the only selection that is taken into account
608 * is egcCompressedX, but other selections should be added when
609 * e.g. writing energies is implemented.
611 static void add_selection_groups(gmx_tng_trajectory_t gmx_tng
,
612 const gmx_mtop_t
*mtop
)
614 const t_atoms
*atoms
;
616 const t_resinfo
*resInfo
;
617 const t_ilist
*ilist
;
620 tng_molecule_t mol
, iterMol
;
627 tng_trajectory_t tng
= gmx_tng
->tng
;
629 /* TODO: When the TNG molecules block is more flexible TNG selection
630 * groups should not need all atoms specified. It should be possible
631 * just to specify what molecules are selected (and/or which atoms in
632 * the molecule) and how many (if applicable). */
634 /* If no atoms are selected we do not need to create a
635 * TNG selection group molecule. */
636 if (mtop
->groups
.ngrpnr
[egcCompressedX
] == 0)
641 /* If all atoms are selected we do not have to create a selection
642 * group molecule in the TNG molecule system. */
643 if (all_atoms_selected(mtop
, egcCompressedX
))
648 /* The name of the TNG molecule containing the selection group is the
649 * same as the name of the selection group. */
650 nameIndex
= *mtop
->groups
.grps
[egcCompressedX
].nm_ind
;
651 groupName
= *mtop
->groups
.grpname
[nameIndex
];
653 tng_molecule_alloc(tng
, &mol
);
654 tng_molecule_name_set(tng
, mol
, groupName
);
655 tng_molecule_chain_add(tng
, mol
, "", &chain
);
657 for (const gmx_molblock_t
&molBlock
: mtop
->molblock
)
659 const gmx_moltype_t
&molType
= mtop
->moltype
[molBlock
.type
];
661 atoms
= &molType
.atoms
;
663 for (int j
= 0; j
< molBlock
.nmol
; j
++)
665 bool bAtomsAdded
= FALSE
;
666 for (int atomIndex
= 0; atomIndex
< atoms
->nr
; atomIndex
++, i
++)
671 if (ggrpnr(&mtop
->groups
, egcCompressedX
, i
) != 0)
675 at
= &atoms
->atom
[atomIndex
];
678 resInfo
= &atoms
->resinfo
[at
->resind
];
679 /* FIXME: When TNG supports both residue index and residue
680 * number the latter should be used. */
681 res_name
= *resInfo
->name
;
682 res_id
= at
->resind
+ 1;
686 res_name
= (char *)"";
689 if (tng_chain_residue_find(tng
, chain
, res_name
, res_id
, &res
)
692 /* Since there is ONE chain for selection groups do not keep the
693 * original residue IDs - otherwise there might be conflicts. */
694 tng_chain_residue_add(tng
, chain
, res_name
, &res
);
696 tng_residue_atom_w_id_add(tng
, res
, *(atoms
->atomname
[atomIndex
]),
697 *(atoms
->atomtype
[atomIndex
]),
698 atom_offset
+ atomIndex
, &atom
);
704 for (int k
= 0; k
< F_NRE
; k
++)
708 ilist
= &molType
.ilist
[k
];
712 while (l
< ilist
->nr
)
715 atom1
= ilist
->iatoms
[l
] + atom_offset
;
716 atom2
= ilist
->iatoms
[l
+1] + atom_offset
;
717 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom1
) == 0 &&
718 ggrpnr(&mtop
->groups
, egcCompressedX
, atom2
) == 0)
720 tng_molecule_bond_add(tng
, mol
, ilist
->iatoms
[l
],
721 ilist
->iatoms
[l
+1], &tngBond
);
728 /* Settle is described using three atoms */
729 ilist
= &molType
.ilist
[F_SETTLE
];
733 while (l
< ilist
->nr
)
735 int atom1
, atom2
, atom3
;
736 atom1
= ilist
->iatoms
[l
] + atom_offset
;
737 atom2
= ilist
->iatoms
[l
+1] + atom_offset
;
738 atom3
= ilist
->iatoms
[l
+2] + atom_offset
;
739 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom1
) == 0)
741 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom2
) == 0)
743 tng_molecule_bond_add(tng
, mol
, atom1
,
746 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom3
) == 0)
748 tng_molecule_bond_add(tng
, mol
, atom1
,
756 atom_offset
+= atoms
->nr
;
759 tng_molecule_existing_add(tng
, &mol
);
760 tng_molecule_cnt_set(tng
, mol
, 1);
761 tng_num_molecule_types_get(tng
, &nMols
);
762 for (gmx_int64_t k
= 0; k
< nMols
; k
++)
764 tng_molecule_of_index_get(tng
, k
, &iterMol
);
769 tng_molecule_cnt_set(tng
, iterMol
, 0);
774 void gmx_tng_set_compression_precision(gmx_tng_trajectory_t gmx_tng
,
778 tng_compression_precision_set(gmx_tng
->tng
, prec
);
780 GMX_UNUSED_VALUE(gmx_tng
);
781 GMX_UNUSED_VALUE(prec
);
785 void gmx_tng_prepare_low_prec_writing(gmx_tng_trajectory_t gmx_tng
,
786 const gmx_mtop_t
*mtop
,
787 const t_inputrec
*ir
)
790 gmx_tng_add_mtop(gmx_tng
, mtop
);
791 add_selection_groups(gmx_tng
, mtop
);
792 set_writing_intervals(gmx_tng
, TRUE
, ir
);
793 tng_time_per_frame_set(gmx_tng
->tng
, ir
->delta_t
* PICO
);
794 gmx_tng
->timePerFrameIsSet
= true;
795 gmx_tng_set_compression_precision(gmx_tng
, ir
->x_compression_precision
);
797 GMX_UNUSED_VALUE(gmx_tng
);
798 GMX_UNUSED_VALUE(mtop
);
799 GMX_UNUSED_VALUE(ir
);
803 void gmx_fwrite_tng(gmx_tng_trajectory_t gmx_tng
,
804 const gmx_bool bUseLossyCompression
,
806 real elapsedPicoSeconds
,
815 typedef tng_function_status (*write_data_func_pointer
)(tng_trajectory_t
,
825 static write_data_func_pointer write_data
= tng_util_generic_with_time_double_write
;
827 static write_data_func_pointer write_data
= tng_util_generic_with_time_write
;
829 double elapsedSeconds
= elapsedPicoSeconds
* PICO
;
830 gmx_int64_t nParticles
;
836 /* This function might get called when the type of the
837 compressed trajectory is actually XTC. So we exit and move
841 tng_trajectory_t tng
= gmx_tng
->tng
;
843 // While the GROMACS interface to this routine specifies 'step', TNG itself
844 // only uses 'frame index' internally, although it suggests that it's a good
845 // idea to let this represent the step rather than just frame numbers.
847 // However, the way the GROMACS interface works, there are cases where
848 // the step is simply not set, so all frames rather have step=0.
849 // If we call the lower-level TNG routines with the same frame index
850 // (which is set from the step) multiple times, things will go very
851 // wrong (overwritten data).
853 // To avoid this, we require the step value to always be larger than the
854 // previous value, and if this is not true we simply set it to a value
855 // one unit larger than the previous step.
857 // Note: We do allow the user to specify any crazy value the want for the
858 // first step. If it is "not set", GROMACS will have used the value 0.
859 if (gmx_tng
->lastStepDataIsValid
&& step
<= gmx_tng
->lastStep
)
861 step
= gmx_tng
->lastStep
+ 1;
864 // Now that we have fixed the potentially incorrect step, we can also set
865 // the time per frame if it was not already done, and we have
866 // step/time information from the last frame so it is possible to calculate it.
867 // Note that we do allow the time to be the same for all steps, or even
868 // decreasing. In that case we will simply say the time per step is 0
869 // or negative. A bit strange, but a correct representation of the data :-)
870 if (!gmx_tng
->timePerFrameIsSet
&& gmx_tng
->lastTimeDataIsValid
&& gmx_tng
->lastStepDataIsValid
)
872 double deltaTime
= elapsedSeconds
- gmx_tng
->lastTime
;
873 std::int64_t deltaStep
= step
- gmx_tng
->lastStep
;
874 tng_time_per_frame_set(tng
, deltaTime
/ deltaStep
);
875 gmx_tng
->timePerFrameIsSet
= true;
878 tng_num_particles_get(tng
, &nParticles
);
879 if (nAtoms
!= (int)nParticles
)
881 tng_implicit_num_particles_set(tng
, nAtoms
);
884 if (bUseLossyCompression
)
886 compression
= TNG_TNG_COMPRESSION
;
890 compression
= TNG_GZIP_COMPRESSION
;
893 /* The writing is done using write_data, which writes float or double
894 * depending on the GROMACS compilation. */
897 GMX_ASSERT(box
, "Need a non-NULL box if positions are written");
899 if (write_data(tng
, step
, elapsedSeconds
,
900 reinterpret_cast<const real
*>(x
),
901 3, TNG_TRAJ_POSITIONS
, "POSITIONS",
902 TNG_PARTICLE_BLOCK_DATA
,
903 compression
) != TNG_SUCCESS
)
905 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
911 if (write_data(tng
, step
, elapsedSeconds
,
912 reinterpret_cast<const real
*>(v
),
913 3, TNG_TRAJ_VELOCITIES
, "VELOCITIES",
914 TNG_PARTICLE_BLOCK_DATA
,
915 compression
) != TNG_SUCCESS
)
917 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
923 /* TNG-MF1 compression only compresses positions and velocities. Use lossless
924 * compression for forces regardless of output mode */
925 if (write_data(tng
, step
, elapsedSeconds
,
926 reinterpret_cast<const real
*>(f
),
927 3, TNG_TRAJ_FORCES
, "FORCES",
928 TNG_PARTICLE_BLOCK_DATA
,
929 TNG_GZIP_COMPRESSION
) != TNG_SUCCESS
)
931 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
937 /* TNG-MF1 compression only compresses positions and velocities. Use lossless
938 * compression for box shape regardless of output mode */
939 if (write_data(tng
, step
, elapsedSeconds
,
940 reinterpret_cast<const real
*>(box
),
941 9, TNG_TRAJ_BOX_SHAPE
, "BOX SHAPE",
942 TNG_NON_PARTICLE_BLOCK_DATA
,
943 TNG_GZIP_COMPRESSION
) != TNG_SUCCESS
)
945 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
951 /* TNG-MF1 compression only compresses positions and velocities. Use lossless
952 * compression for lambda regardless of output mode */
953 if (write_data(tng
, step
, elapsedSeconds
,
954 reinterpret_cast<const real
*>(&lambda
),
955 1, TNG_GMX_LAMBDA
, "LAMBDAS",
956 TNG_NON_PARTICLE_BLOCK_DATA
,
957 TNG_GZIP_COMPRESSION
) != TNG_SUCCESS
)
959 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
963 // Update the data in the wrapper
965 gmx_tng
->lastStepDataIsValid
= true;
966 gmx_tng
->lastStep
= step
;
967 gmx_tng
->lastTimeDataIsValid
= true;
968 gmx_tng
->lastTime
= elapsedSeconds
;
970 GMX_UNUSED_VALUE(gmx_tng
);
971 GMX_UNUSED_VALUE(bUseLossyCompression
);
972 GMX_UNUSED_VALUE(step
);
973 GMX_UNUSED_VALUE(elapsedPicoSeconds
);
974 GMX_UNUSED_VALUE(lambda
);
975 GMX_UNUSED_VALUE(box
);
976 GMX_UNUSED_VALUE(nAtoms
);
983 void fflush_tng(gmx_tng_trajectory_t gmx_tng
)
990 tng_frame_set_premature_write(gmx_tng
->tng
, TNG_USE_HASH
);
992 GMX_UNUSED_VALUE(gmx_tng
);
996 float gmx_tng_get_time_of_final_frame(gmx_tng_trajectory_t gmx_tng
)
1002 tng_trajectory_t tng
= gmx_tng
->tng
;
1004 tng_num_frames_get(tng
, &nFrames
);
1005 tng_util_time_of_frame_get(tng
, nFrames
- 1, &time
);
1007 fTime
= time
/ PICO
;
1010 GMX_UNUSED_VALUE(gmx_tng
);
1015 void gmx_prepare_tng_writing(const char *filename
,
1017 gmx_tng_trajectory_t
*gmx_tng_input
,
1018 gmx_tng_trajectory_t
*gmx_tng_output
,
1020 const gmx_mtop_t
*mtop
,
1022 const char *indexGroupName
)
1025 tng_trajectory_t
*input
= (gmx_tng_input
&& *gmx_tng_input
) ? &(*gmx_tng_input
)->tng
: nullptr;
1026 /* FIXME after 5.0: Currently only standard block types are read */
1027 const int defaultNumIds
= 5;
1028 static gmx_int64_t fallbackIds
[defaultNumIds
] =
1030 TNG_TRAJ_BOX_SHAPE
, TNG_TRAJ_POSITIONS
,
1031 TNG_TRAJ_VELOCITIES
, TNG_TRAJ_FORCES
,
1034 static char fallbackNames
[defaultNumIds
][32] =
1036 "BOX SHAPE", "POSITIONS", "VELOCITIES",
1040 typedef tng_function_status (*set_writing_interval_func_pointer
)(tng_trajectory_t
,
1048 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_double_set
;
1050 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_set
;
1053 gmx_tng_open(filename
, mode
, gmx_tng_output
);
1054 tng_trajectory_t
*output
= &(*gmx_tng_output
)->tng
;
1056 /* Do we have an input file in TNG format? If so, then there's
1057 more data we can copy over, rather than having to improvise. */
1058 if (gmx_tng_input
&& *gmx_tng_input
)
1060 /* Set parameters (compression, time per frame, molecule
1061 * information, number of frames per frame set and writing
1062 * intervals of positions, box shape and lambdas) of the
1063 * output tng container based on their respective values int
1064 * the input tng container */
1065 double time
, compression_precision
;
1066 gmx_int64_t n_frames_per_frame_set
, interval
= -1;
1068 tng_compression_precision_get(*input
, &compression_precision
);
1069 tng_compression_precision_set(*output
, compression_precision
);
1070 // TODO make this configurable in a future version
1071 char compression_type
= TNG_TNG_COMPRESSION
;
1073 tng_molecule_system_copy(*input
, *output
);
1075 tng_time_per_frame_get(*input
, &time
);
1076 tng_time_per_frame_set(*output
, time
);
1077 // Since we have copied the value from the input TNG we should not change it again
1078 (*gmx_tng_output
)->timePerFrameIsSet
= true;
1080 tng_num_frames_per_frame_set_get(*input
, &n_frames_per_frame_set
);
1081 tng_num_frames_per_frame_set_set(*output
, n_frames_per_frame_set
);
1083 for (int i
= 0; i
< defaultNumIds
; i
++)
1085 if (tng_data_get_stride_length(*input
, fallbackIds
[i
], -1, &interval
)
1088 switch (fallbackIds
[i
])
1090 case TNG_TRAJ_POSITIONS
:
1091 case TNG_TRAJ_VELOCITIES
:
1092 set_writing_interval(*output
, interval
, 3, fallbackIds
[i
],
1093 fallbackNames
[i
], TNG_PARTICLE_BLOCK_DATA
,
1096 case TNG_TRAJ_FORCES
:
1097 set_writing_interval(*output
, interval
, 3, fallbackIds
[i
],
1098 fallbackNames
[i
], TNG_PARTICLE_BLOCK_DATA
,
1099 TNG_GZIP_COMPRESSION
);
1101 case TNG_TRAJ_BOX_SHAPE
:
1102 set_writing_interval(*output
, interval
, 9, fallbackIds
[i
],
1103 fallbackNames
[i
], TNG_NON_PARTICLE_BLOCK_DATA
,
1104 TNG_GZIP_COMPRESSION
);
1105 (*gmx_tng_output
)->boxOutputInterval
= interval
;
1107 case TNG_GMX_LAMBDA
:
1108 set_writing_interval(*output
, interval
, 1, fallbackIds
[i
],
1109 fallbackNames
[i
], TNG_NON_PARTICLE_BLOCK_DATA
,
1110 TNG_GZIP_COMPRESSION
);
1111 (*gmx_tng_output
)->lambdaOutputInterval
= interval
;
1122 /* TODO after trjconv is modularized: fix this so the user can
1123 change precision when they are doing an operation where
1124 this makes sense, and not otherwise.
1126 char compression = bUseLossyCompression ? TNG_TNG_COMPRESSION : TNG_GZIP_COMPRESSION;
1127 gmx_tng_set_compression_precision(*output, ndec2prec(nDecimalsOfPrecision));
1129 gmx_tng_add_mtop(*gmx_tng_output
, mtop
);
1130 tng_num_frames_per_frame_set_set(*output
, 1);
1133 if (index
&& nAtoms
> 0)
1135 gmx_tng_setup_atom_subgroup(*gmx_tng_output
, nAtoms
, index
, indexGroupName
);
1138 /* If for some reason there are more requested atoms than there are atoms in the
1139 * molecular system create a number of implicit atoms (without atom data) to
1140 * compensate for that. */
1143 tng_implicit_num_particles_set(*output
, nAtoms
);
1146 GMX_UNUSED_VALUE(filename
);
1147 GMX_UNUSED_VALUE(mode
);
1148 GMX_UNUSED_VALUE(gmx_tng_input
);
1149 GMX_UNUSED_VALUE(gmx_tng_output
);
1150 GMX_UNUSED_VALUE(nAtoms
);
1151 GMX_UNUSED_VALUE(mtop
);
1152 GMX_UNUSED_VALUE(index
);
1153 GMX_UNUSED_VALUE(indexGroupName
);
1157 void gmx_write_tng_from_trxframe(gmx_tng_trajectory_t gmx_tng_output
,
1158 const t_trxframe
*frame
,
1164 natoms
= frame
->natoms
;
1166 gmx_fwrite_tng(gmx_tng_output
,
1177 GMX_UNUSED_VALUE(gmx_tng_output
);
1178 GMX_UNUSED_VALUE(frame
);
1179 GMX_UNUSED_VALUE(natoms
);
1188 convert_array_to_real_array(void *from
,
1193 const char datatype
)
1197 const bool useDouble
= GMX_DOUBLE
;
1200 case TNG_FLOAT_DATA
:
1205 memcpy(to
, from
, nValues
* sizeof(real
) * nAtoms
);
1209 for (i
= 0; i
< nAtoms
; i
++)
1211 for (j
= 0; j
< nValues
; j
++)
1213 to
[i
*nValues
+j
] = reinterpret_cast<float *>(from
)[i
*nValues
+j
] * fact
;
1220 for (i
= 0; i
< nAtoms
; i
++)
1222 for (j
= 0; j
< nValues
; j
++)
1224 to
[i
*nValues
+j
] = reinterpret_cast<float *>(from
)[i
*nValues
+j
] * fact
;
1230 for (i
= 0; i
< nAtoms
; i
++)
1232 for (j
= 0; j
< nValues
; j
++)
1234 to
[i
*nValues
+j
] = reinterpret_cast<gmx_int64_t
*>(from
)[i
*nValues
+j
] * fact
;
1238 case TNG_DOUBLE_DATA
:
1239 if (sizeof(real
) == sizeof(double))
1243 memcpy(to
, from
, nValues
* sizeof(real
) * nAtoms
);
1247 for (i
= 0; i
< nAtoms
; i
++)
1249 for (j
= 0; j
< nValues
; j
++)
1251 to
[i
*nValues
+j
] = reinterpret_cast<double *>(from
)[i
*nValues
+j
] * fact
;
1258 for (i
= 0; i
< nAtoms
; i
++)
1260 for (j
= 0; j
< nValues
; j
++)
1262 to
[i
*nValues
+j
] = reinterpret_cast<double *>(from
)[i
*nValues
+j
] * fact
;
1268 gmx_incons("Illegal datatype when converting values to a real array!");
1274 real
getDistanceScaleFactor(gmx_tng_trajectory_t in
)
1276 gmx_int64_t exp
= -1;
1277 real distanceScaleFactor
;
1279 // TODO Hopefully, TNG 2.0 will do this kind of thing for us
1280 tng_distance_unit_exponential_get(in
->tng
, &exp
);
1282 // GROMACS expects distances in nm
1286 distanceScaleFactor
= NANO
/NANO
;
1289 distanceScaleFactor
= NANO
/ANGSTROM
;
1292 distanceScaleFactor
= pow(10.0, exp
+ 9.0);
1295 return distanceScaleFactor
;
1301 void gmx_tng_setup_atom_subgroup(gmx_tng_trajectory_t gmx_tng
,
1307 gmx_int64_t nAtoms
, cnt
, nMols
;
1308 tng_molecule_t mol
, iterMol
;
1312 tng_function_status stat
;
1313 tng_trajectory_t tng
= gmx_tng
->tng
;
1315 tng_num_particles_get(tng
, &nAtoms
);
1322 stat
= tng_molecule_find(tng
, name
, -1, &mol
);
1323 if (stat
== TNG_SUCCESS
)
1325 tng_molecule_num_atoms_get(tng
, mol
, &nAtoms
);
1326 tng_molecule_cnt_get(tng
, mol
, &cnt
);
1336 if (stat
== TNG_FAILURE
)
1338 /* The indexed atoms are added to one separate molecule. */
1339 tng_molecule_alloc(tng
, &mol
);
1340 tng_molecule_name_set(tng
, mol
, name
);
1341 tng_molecule_chain_add(tng
, mol
, "", &chain
);
1343 for (int i
= 0; i
< nind
; i
++)
1345 char temp_name
[256], temp_type
[256];
1347 /* Try to retrieve the residue name of the atom */
1348 stat
= tng_residue_name_of_particle_nr_get(tng
, ind
[i
], temp_name
, 256);
1349 if (stat
!= TNG_SUCCESS
)
1351 temp_name
[0] = '\0';
1353 /* Check if the molecule of the selection already contains this residue */
1354 if (tng_chain_residue_find(tng
, chain
, temp_name
, -1, &res
)
1357 tng_chain_residue_add(tng
, chain
, temp_name
, &res
);
1359 /* Try to find the original name and type of the atom */
1360 stat
= tng_atom_name_of_particle_nr_get(tng
, ind
[i
], temp_name
, 256);
1361 if (stat
!= TNG_SUCCESS
)
1363 temp_name
[0] = '\0';
1365 stat
= tng_atom_type_of_particle_nr_get(tng
, ind
[i
], temp_type
, 256);
1366 if (stat
!= TNG_SUCCESS
)
1368 temp_type
[0] = '\0';
1370 tng_residue_atom_w_id_add(tng
, res
, temp_name
, temp_type
, ind
[i
], &atom
);
1372 tng_molecule_existing_add(tng
, &mol
);
1374 /* Set the count of the molecule containing the selected atoms to 1 and all
1375 * other molecules to 0 */
1376 tng_molecule_cnt_set(tng
, mol
, 1);
1377 tng_num_molecule_types_get(tng
, &nMols
);
1378 for (gmx_int64_t k
= 0; k
< nMols
; k
++)
1380 tng_molecule_of_index_get(tng
, k
, &iterMol
);
1385 tng_molecule_cnt_set(tng
, iterMol
, 0);
1388 GMX_UNUSED_VALUE(gmx_tng
);
1389 GMX_UNUSED_VALUE(nind
);
1390 GMX_UNUSED_VALUE(ind
);
1391 GMX_UNUSED_VALUE(name
);
1395 /* TODO: If/when TNG acquires the ability to copy data blocks without
1396 * uncompressing them, then this implemenation should be reconsidered.
1397 * Ideally, gmx trjconv -f a.tng -o b.tng -b 10 -e 20 would be fast
1398 * and lose no information. */
1399 gmx_bool
gmx_read_next_tng_frame(gmx_tng_trajectory_t gmx_tng_input
,
1401 gmx_int64_t
*requestedIds
,
1402 int numRequestedIds
)
1405 tng_trajectory_t input
= gmx_tng_input
->tng
;
1406 gmx_bool bOK
= TRUE
;
1407 tng_function_status stat
;
1408 gmx_int64_t numberOfAtoms
= -1, frameNumber
= -1;
1409 gmx_int64_t nBlocks
, blockId
, *blockIds
= nullptr, codecId
;
1411 void *values
= nullptr;
1412 double frameTime
= -1.0;
1413 int size
, blockDependency
;
1415 const int defaultNumIds
= 5;
1416 static gmx_int64_t fallbackRequestedIds
[defaultNumIds
] =
1418 TNG_TRAJ_BOX_SHAPE
, TNG_TRAJ_POSITIONS
,
1419 TNG_TRAJ_VELOCITIES
, TNG_TRAJ_FORCES
,
1426 fr
->bLambda
= FALSE
;
1434 /* If no specific IDs were requested read all block types that can
1435 * currently be interpreted */
1436 if (!requestedIds
|| numRequestedIds
== 0)
1438 numRequestedIds
= defaultNumIds
;
1439 requestedIds
= fallbackRequestedIds
;
1442 stat
= tng_num_particles_get(input
, &numberOfAtoms
);
1443 if (stat
!= TNG_SUCCESS
)
1445 gmx_file("Cannot determine number of atoms from TNG file.");
1447 fr
->natoms
= numberOfAtoms
;
1449 bool nextFrameExists
= gmx_get_tng_data_block_types_of_next_frame(gmx_tng_input
,
1456 gmx::unique_cptr
<gmx_int64_t
, gmx::free_wrapper
> blockIdsGuard(blockIds
);
1457 if (!nextFrameExists
)
1467 for (gmx_int64_t i
= 0; i
< nBlocks
; i
++)
1469 blockId
= blockIds
[i
];
1470 tng_data_block_dependency_get(input
, blockId
, &blockDependency
);
1471 if (blockDependency
& TNG_PARTICLE_DEPENDENT
)
1473 stat
= tng_util_particle_data_next_frame_read(input
,
1482 stat
= tng_util_non_particle_data_next_frame_read(input
,
1489 if (stat
== TNG_CRITICAL
)
1491 gmx_file("Cannot read positions from TNG file.");
1494 else if (stat
== TNG_FAILURE
)
1500 case TNG_TRAJ_BOX_SHAPE
:
1504 size
= sizeof(gmx_int64_t
);
1506 case TNG_FLOAT_DATA
:
1507 size
= sizeof(float);
1509 case TNG_DOUBLE_DATA
:
1510 size
= sizeof(double);
1513 gmx_incons("Illegal datatype of box shape values!");
1515 for (int i
= 0; i
< DIM
; i
++)
1517 convert_array_to_real_array(reinterpret_cast<char *>(values
) + size
* i
* DIM
,
1518 reinterpret_cast<real
*>(fr
->box
[i
]),
1519 getDistanceScaleFactor(gmx_tng_input
),
1526 case TNG_TRAJ_POSITIONS
:
1527 srenew(fr
->x
, fr
->natoms
);
1528 convert_array_to_real_array(values
,
1529 reinterpret_cast<real
*>(fr
->x
),
1530 getDistanceScaleFactor(gmx_tng_input
),
1535 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &prec
);
1536 /* This must be updated if/when more lossy compression methods are added */
1537 if (codecId
== TNG_TNG_COMPRESSION
)
1543 case TNG_TRAJ_VELOCITIES
:
1544 srenew(fr
->v
, fr
->natoms
);
1545 convert_array_to_real_array(values
,
1547 getDistanceScaleFactor(gmx_tng_input
),
1552 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &prec
);
1553 /* This must be updated if/when more lossy compression methods are added */
1554 if (codecId
== TNG_TNG_COMPRESSION
)
1560 case TNG_TRAJ_FORCES
:
1561 srenew(fr
->f
, fr
->natoms
);
1562 convert_array_to_real_array(values
,
1563 reinterpret_cast<real
*>(fr
->f
),
1564 getDistanceScaleFactor(gmx_tng_input
),
1570 case TNG_GMX_LAMBDA
:
1573 case TNG_FLOAT_DATA
:
1574 fr
->lambda
= *(reinterpret_cast<float *>(values
));
1576 case TNG_DOUBLE_DATA
:
1577 fr
->lambda
= *(reinterpret_cast<double *>(values
));
1580 gmx_incons("Illegal datatype lambda value!");
1585 gmx_warning("Illegal block type! Currently GROMACS tools can only handle certain data types. Skipping block.");
1587 /* values does not have to be freed before reading next frame. It will
1588 * be reallocated if it is not NULL. */
1591 fr
->step
= frameNumber
;
1594 // Convert the time to ps
1595 fr
->time
= frameTime
/ PICO
;
1596 fr
->bTime
= (frameTime
> 0);
1598 // TODO This does not leak, but is not exception safe.
1599 /* values must be freed before leaving this function */
1604 GMX_UNUSED_VALUE(gmx_tng_input
);
1605 GMX_UNUSED_VALUE(fr
);
1606 GMX_UNUSED_VALUE(requestedIds
);
1607 GMX_UNUSED_VALUE(numRequestedIds
);
1612 void gmx_print_tng_molecule_system(gmx_tng_trajectory_t gmx_tng_input
,
1616 gmx_int64_t nMolecules
, nChains
, nResidues
, nAtoms
, nFramesRead
;
1617 gmx_int64_t strideLength
, nParticlesRead
, nValuesPerFrameRead
, *molCntList
;
1618 tng_molecule_t molecule
;
1620 tng_residue_t residue
;
1622 tng_function_status stat
;
1626 void *data
= nullptr;
1627 std::vector
<real
> atomCharges
;
1628 std::vector
<real
> atomMasses
;
1629 tng_trajectory_t input
= gmx_tng_input
->tng
;
1631 tng_num_molecule_types_get(input
, &nMolecules
);
1632 tng_molecule_cnt_list_get(input
, &molCntList
);
1633 /* Can the number of particles change in the trajectory or is it constant? */
1634 tng_num_particles_variable_get(input
, &varNAtoms
);
1636 for (gmx_int64_t i
= 0; i
< nMolecules
; i
++)
1638 tng_molecule_of_index_get(input
, i
, &molecule
);
1639 tng_molecule_name_get(input
, molecule
, str
, 256);
1640 if (varNAtoms
== TNG_CONSTANT_N_ATOMS
)
1642 if ((int)molCntList
[i
] == 0)
1646 fprintf(stream
, "Molecule: %s, count: %d\n", str
, (int)molCntList
[i
]);
1650 fprintf(stream
, "Molecule: %s\n", str
);
1652 tng_molecule_num_chains_get(input
, molecule
, &nChains
);
1655 for (gmx_int64_t j
= 0; j
< nChains
; j
++)
1657 tng_molecule_chain_of_index_get(input
, molecule
, j
, &chain
);
1658 tng_chain_name_get(input
, chain
, str
, 256);
1659 fprintf(stream
, "\tChain: %s\n", str
);
1660 tng_chain_num_residues_get(input
, chain
, &nResidues
);
1661 for (gmx_int64_t k
= 0; k
< nResidues
; k
++)
1663 tng_chain_residue_of_index_get(input
, chain
, k
, &residue
);
1664 tng_residue_name_get(input
, residue
, str
, 256);
1665 fprintf(stream
, "\t\tResidue: %s\n", str
);
1666 tng_residue_num_atoms_get(input
, residue
, &nAtoms
);
1667 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
1669 tng_residue_atom_of_index_get(input
, residue
, l
, &atom
);
1670 tng_atom_name_get(input
, atom
, str
, 256);
1671 fprintf(stream
, "\t\t\tAtom: %s", str
);
1672 tng_atom_type_get(input
, atom
, str
, 256);
1673 fprintf(stream
, " (%s)\n", str
);
1678 /* It is possible to have a molecule without chains, in which case
1679 * residues in the molecule can be iterated through without going
1680 * through chains. */
1683 tng_molecule_num_residues_get(input
, molecule
, &nResidues
);
1686 for (gmx_int64_t k
= 0; k
< nResidues
; k
++)
1688 tng_molecule_residue_of_index_get(input
, molecule
, k
, &residue
);
1689 tng_residue_name_get(input
, residue
, str
, 256);
1690 fprintf(stream
, "\t\tResidue: %s\n", str
);
1691 tng_residue_num_atoms_get(input
, residue
, &nAtoms
);
1692 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
1694 tng_residue_atom_of_index_get(input
, residue
, l
, &atom
);
1695 tng_atom_name_get(input
, atom
, str
, 256);
1696 fprintf(stream
, "\t\t\tAtom: %s", str
);
1697 tng_atom_type_get(input
, atom
, str
, 256);
1698 fprintf(stream
, " (%s)\n", str
);
1704 tng_molecule_num_atoms_get(input
, molecule
, &nAtoms
);
1705 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
1707 tng_molecule_atom_of_index_get(input
, molecule
, l
, &atom
);
1708 tng_atom_name_get(input
, atom
, str
, 256);
1709 fprintf(stream
, "\t\t\tAtom: %s", str
);
1710 tng_atom_type_get(input
, atom
, str
, 256);
1711 fprintf(stream
, " (%s)\n", str
);
1717 tng_num_particles_get(input
, &nAtoms
);
1718 stat
= tng_particle_data_vector_get(input
, TNG_TRAJ_PARTIAL_CHARGES
, &data
, &nFramesRead
,
1719 &strideLength
, &nParticlesRead
,
1720 &nValuesPerFrameRead
, &datatype
);
1721 if (stat
== TNG_SUCCESS
)
1723 atomCharges
.resize(nAtoms
);
1724 convert_array_to_real_array(data
,
1731 fprintf(stream
, "Atom Charges (%d):\n", int(nAtoms
));
1732 for (gmx_int64_t i
= 0; i
< nAtoms
; i
+= 10)
1734 fprintf(stream
, "Atom Charges [%8d-]=[", int(i
));
1735 for (gmx_int64_t j
= 0; (j
< 10 && i
+ j
< nAtoms
); j
++)
1737 fprintf(stream
, " %12.5e", atomCharges
[i
+ j
]);
1739 fprintf(stream
, "]\n");
1743 stat
= tng_particle_data_vector_get(input
, TNG_TRAJ_MASSES
, &data
, &nFramesRead
,
1744 &strideLength
, &nParticlesRead
,
1745 &nValuesPerFrameRead
, &datatype
);
1746 if (stat
== TNG_SUCCESS
)
1748 atomMasses
.resize(nAtoms
);
1749 convert_array_to_real_array(data
,
1756 fprintf(stream
, "Atom Masses (%d):\n", int(nAtoms
));
1757 for (gmx_int64_t i
= 0; i
< nAtoms
; i
+= 10)
1759 fprintf(stream
, "Atom Masses [%8d-]=[", int(i
));
1760 for (gmx_int64_t j
= 0; (j
< 10 && i
+ j
< nAtoms
); j
++)
1762 fprintf(stream
, " %12.5e", atomMasses
[i
+ j
]);
1764 fprintf(stream
, "]\n");
1770 GMX_UNUSED_VALUE(gmx_tng_input
);
1771 GMX_UNUSED_VALUE(stream
);
1775 gmx_bool
gmx_get_tng_data_block_types_of_next_frame(gmx_tng_trajectory_t gmx_tng_input
,
1778 gmx_int64_t
*requestedIds
,
1779 gmx_int64_t
*nextFrame
,
1780 gmx_int64_t
*nBlocks
,
1781 gmx_int64_t
**blockIds
)
1784 tng_function_status stat
;
1785 tng_trajectory_t input
= gmx_tng_input
->tng
;
1787 stat
= tng_util_trajectory_next_frame_present_data_blocks_find(input
, frame
,
1788 nRequestedIds
, requestedIds
,
1792 if (stat
== TNG_CRITICAL
)
1794 gmx_file("Cannot read TNG file. Cannot find data blocks of next frame.");
1796 else if (stat
== TNG_FAILURE
)
1802 GMX_UNUSED_VALUE(gmx_tng_input
);
1803 GMX_UNUSED_VALUE(frame
);
1804 GMX_UNUSED_VALUE(nRequestedIds
);
1805 GMX_UNUSED_VALUE(requestedIds
);
1806 GMX_UNUSED_VALUE(nextFrame
);
1807 GMX_UNUSED_VALUE(nBlocks
);
1808 GMX_UNUSED_VALUE(blockIds
);
1813 gmx_bool
gmx_get_tng_data_next_frame_of_block_type(gmx_tng_trajectory_t gmx_tng_input
,
1814 gmx_int64_t blockId
,
1816 gmx_int64_t
*frameNumber
,
1818 gmx_int64_t
*nValuesPerFrame
,
1819 gmx_int64_t
*nAtoms
,
1826 tng_function_status stat
;
1828 gmx_int64_t codecId
;
1829 int blockDependency
;
1830 void *data
= nullptr;
1832 tng_trajectory_t input
= gmx_tng_input
->tng
;
1834 stat
= tng_data_block_name_get(input
, blockId
, name
, maxLen
);
1835 if (stat
!= TNG_SUCCESS
)
1837 gmx_file("Cannot read next frame of TNG file");
1839 stat
= tng_data_block_dependency_get(input
, blockId
, &blockDependency
);
1840 if (stat
!= TNG_SUCCESS
)
1842 gmx_file("Cannot read next frame of TNG file");
1844 if (blockDependency
& TNG_PARTICLE_DEPENDENT
)
1846 tng_num_particles_get(input
, nAtoms
);
1847 stat
= tng_util_particle_data_next_frame_read(input
,
1856 *nAtoms
= 1; /* There are not actually any atoms, but it is used for
1857 allocating memory */
1858 stat
= tng_util_non_particle_data_next_frame_read(input
,
1865 if (stat
== TNG_CRITICAL
)
1867 gmx_file("Cannot read next frame of TNG file");
1869 if (stat
== TNG_FAILURE
)
1875 stat
= tng_data_block_num_values_per_frame_get(input
, blockId
, nValuesPerFrame
);
1876 if (stat
!= TNG_SUCCESS
)
1878 gmx_file("Cannot read next frame of TNG file");
1880 srenew(*values
, sizeof(real
) * *nValuesPerFrame
* *nAtoms
);
1881 convert_array_to_real_array(data
,
1883 getDistanceScaleFactor(gmx_tng_input
),
1888 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &localPrec
);
1890 /* This must be updated if/when more lossy compression methods are added */
1891 if (codecId
!= TNG_TNG_COMPRESSION
)
1905 GMX_UNUSED_VALUE(gmx_tng_input
);
1906 GMX_UNUSED_VALUE(blockId
);
1907 GMX_UNUSED_VALUE(values
);
1908 GMX_UNUSED_VALUE(frameNumber
);
1909 GMX_UNUSED_VALUE(frameTime
);
1910 GMX_UNUSED_VALUE(nValuesPerFrame
);
1911 GMX_UNUSED_VALUE(nAtoms
);
1912 GMX_UNUSED_VALUE(prec
);
1913 GMX_UNUSED_VALUE(name
);
1914 GMX_UNUSED_VALUE(maxLen
);
1915 GMX_UNUSED_VALUE(bOK
);
1920 int gmx_tng_get_box_output_interval(gmx_tng_trajectory_t gmx_tng
)
1923 return gmx_tng
->boxOutputInterval
;
1925 GMX_UNUSED_VALUE(gmx_tng
);
1929 int gmx_tng_get_lambda_output_interval(gmx_tng_trajectory_t gmx_tng
)
1932 return gmx_tng
->lambdaOutputInterval
;
1934 GMX_UNUSED_VALUE(gmx_tng
);