2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014,2015, 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.
37 #include "tngio_for_tools.h"
44 #include "tng/tng_io.h"
47 #include "gromacs/fileio/tngio.h"
48 #include "gromacs/math/units.h"
49 #include "gromacs/trajectory/trajectoryframe.h"
50 #include "gromacs/utility/basedefinitions.h"
51 #include "gromacs/utility/fatalerror.h"
52 #include "gromacs/utility/smalloc.h"
54 void gmx_prepare_tng_writing(const char *filename
,
56 tng_trajectory_t
*input
,
57 tng_trajectory_t
*output
,
59 const gmx_mtop_t
*mtop
,
61 const char *indexGroupName
)
64 /* FIXME after 5.0: Currently only standard block types are read */
65 const int defaultNumIds
= 5;
66 static gmx_int64_t fallbackIds
[defaultNumIds
] =
68 TNG_TRAJ_BOX_SHAPE
, TNG_TRAJ_POSITIONS
,
69 TNG_TRAJ_VELOCITIES
, TNG_TRAJ_FORCES
,
72 static char fallbackNames
[defaultNumIds
][32] =
74 "BOX SHAPE", "POSITIONS", "VELOCITIES",
78 typedef tng_function_status (*set_writing_interval_func_pointer
)(tng_trajectory_t
,
86 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_double_set
;
88 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_set
;
91 gmx_tng_open(filename
, mode
, output
);
93 /* Do we have an input file in TNG format? If so, then there's
94 more data we can copy over, rather than having to improvise. */
97 /* Set parameters (compression, time per frame, molecule
98 * information, number of frames per frame set and writing
99 * intervals of positions, box shape and lambdas) of the
100 * output tng container based on their respective values int
101 * the input tng container */
102 double time
, compression_precision
;
103 gmx_int64_t n_frames_per_frame_set
, interval
= -1;
105 tng_compression_precision_get(*input
, &compression_precision
);
106 tng_compression_precision_set(*output
, compression_precision
);
107 // TODO make this configurable in a future version
108 char compression_type
= TNG_TNG_COMPRESSION
;
110 tng_molecule_system_copy(*input
, *output
);
112 tng_time_per_frame_get(*input
, &time
);
113 tng_time_per_frame_set(*output
, time
);
115 tng_num_frames_per_frame_set_get(*input
, &n_frames_per_frame_set
);
116 tng_num_frames_per_frame_set_set(*output
, n_frames_per_frame_set
);
118 for (int i
= 0; i
< defaultNumIds
; i
++)
120 if (tng_data_get_stride_length(*input
, fallbackIds
[i
], -1, &interval
)
123 switch (fallbackIds
[i
])
125 case TNG_TRAJ_POSITIONS
:
126 case TNG_TRAJ_VELOCITIES
:
127 set_writing_interval(*output
, interval
, 3, fallbackIds
[i
],
128 fallbackNames
[i
], TNG_PARTICLE_BLOCK_DATA
,
131 case TNG_TRAJ_FORCES
:
132 set_writing_interval(*output
, interval
, 3, fallbackIds
[i
],
133 fallbackNames
[i
], TNG_PARTICLE_BLOCK_DATA
,
134 TNG_GZIP_COMPRESSION
);
136 case TNG_TRAJ_BOX_SHAPE
:
137 set_writing_interval(*output
, interval
, 9, fallbackIds
[i
],
138 fallbackNames
[i
], TNG_NON_PARTICLE_BLOCK_DATA
,
139 TNG_GZIP_COMPRESSION
);
142 set_writing_interval(*output
, interval
, 1, fallbackIds
[i
],
143 fallbackNames
[i
], TNG_NON_PARTICLE_BLOCK_DATA
,
144 TNG_GZIP_COMPRESSION
);
154 /* TODO after trjconv is modularized: fix this so the user can
155 change precision when they are doing an operation where
156 this makes sense, and not otherwise.
158 char compression = bUseLossyCompression ? TNG_TNG_COMPRESSION : TNG_GZIP_COMPRESSION;
159 gmx_tng_set_compression_precision(*output, ndec2prec(nDecimalsOfPrecision));
161 gmx_tng_add_mtop(*output
, mtop
);
162 tng_num_frames_per_frame_set_set(*output
, 1);
165 if (index
&& nAtoms
> 0)
167 gmx_tng_setup_atom_subgroup(*output
, nAtoms
, index
, indexGroupName
);
170 /* If for some reason there are more requested atoms than there are atoms in the
171 * molecular system create a number of implicit atoms (without atom data) to
172 * compensate for that. */
175 tng_implicit_num_particles_set(*output
, nAtoms
);
178 GMX_UNUSED_VALUE(filename
);
179 GMX_UNUSED_VALUE(mode
);
180 GMX_UNUSED_VALUE(input
);
181 GMX_UNUSED_VALUE(output
);
182 GMX_UNUSED_VALUE(nAtoms
);
183 GMX_UNUSED_VALUE(mtop
);
184 GMX_UNUSED_VALUE(index
);
185 GMX_UNUSED_VALUE(indexGroupName
);
189 void gmx_write_tng_from_trxframe(tng_trajectory_t output
,
196 double timePerFrame
= frame
->time
* PICO
/ frame
->step
;
197 tng_time_per_frame_set(output
, timePerFrame
);
201 natoms
= frame
->natoms
;
203 gmx_fwrite_tng(output
,
214 GMX_UNUSED_VALUE(output
);
215 GMX_UNUSED_VALUE(frame
);
216 GMX_UNUSED_VALUE(natoms
);
222 convert_array_to_real_array(void *from
,
234 if (sizeof(real
) == sizeof(float))
238 memcpy(to
, from
, nValues
* sizeof(real
) * nAtoms
);
242 for (i
= 0; i
< nAtoms
; i
++)
244 for (j
= 0; j
< nValues
; j
++)
246 to
[i
*nValues
+j
] = reinterpret_cast<float *>(from
)[i
*nValues
+j
] * fact
;
253 for (i
= 0; i
< nAtoms
; i
++)
255 for (j
= 0; j
< nValues
; j
++)
257 to
[i
*nValues
+j
] = reinterpret_cast<float *>(from
)[i
*nValues
+j
] * fact
;
263 for (i
= 0; i
< nAtoms
; i
++)
265 for (j
= 0; j
< nValues
; j
++)
267 to
[i
*nValues
+j
] = reinterpret_cast<gmx_int64_t
*>(from
)[i
*nValues
+j
] * fact
;
271 case TNG_DOUBLE_DATA
:
272 if (sizeof(real
) == sizeof(double))
276 memcpy(to
, from
, nValues
* sizeof(real
) * nAtoms
);
280 for (i
= 0; i
< nAtoms
; i
++)
282 for (j
= 0; j
< nValues
; j
++)
284 to
[i
*nValues
+j
] = reinterpret_cast<double *>(from
)[i
*nValues
+j
] * fact
;
291 for (i
= 0; i
< nAtoms
; i
++)
293 for (j
= 0; j
< nValues
; j
++)
295 to
[i
*nValues
+j
] = reinterpret_cast<double *>(from
)[i
*nValues
+j
] * fact
;
301 gmx_incons("Illegal datatype when converting values to a real array!");
307 static real
getDistanceScaleFactor(tng_trajectory_t in
)
309 gmx_int64_t exp
= -1;
310 real distanceScaleFactor
;
312 // TODO Hopefully, TNG 2.0 will do this kind of thing for us
313 tng_distance_unit_exponential_get(in
, &exp
);
315 // GROMACS expects distances in nm
319 distanceScaleFactor
= NANO
/NANO
;
322 distanceScaleFactor
= NANO
/ANGSTROM
;
325 distanceScaleFactor
= pow(10.0, exp
+ 9.0);
328 return distanceScaleFactor
;
332 void gmx_tng_setup_atom_subgroup(tng_trajectory_t tng
,
338 gmx_int64_t nAtoms
, cnt
, nMols
;
339 tng_molecule_t mol
, iterMol
;
343 tng_function_status stat
;
345 tng_num_particles_get(tng
, &nAtoms
);
352 stat
= tng_molecule_find(tng
, name
, -1, &mol
);
353 if (stat
== TNG_SUCCESS
)
355 tng_molecule_num_atoms_get(tng
, mol
, &nAtoms
);
356 tng_molecule_cnt_get(tng
, mol
, &cnt
);
366 if (stat
== TNG_FAILURE
)
368 /* The indexed atoms are added to one separate molecule. */
369 tng_molecule_alloc(tng
, &mol
);
370 tng_molecule_name_set(tng
, mol
, name
);
371 tng_molecule_chain_add(tng
, mol
, "", &chain
);
373 for (int i
= 0; i
< nind
; i
++)
375 char temp_name
[256], temp_type
[256];
377 /* Try to retrieve the residue name of the atom */
378 stat
= tng_residue_name_of_particle_nr_get(tng
, ind
[i
], temp_name
, 256);
379 if (stat
!= TNG_SUCCESS
)
383 /* Check if the molecule of the selection already contains this residue */
384 if (tng_chain_residue_find(tng
, chain
, temp_name
, -1, &res
)
387 tng_chain_residue_add(tng
, chain
, temp_name
, &res
);
389 /* Try to find the original name and type of the atom */
390 stat
= tng_atom_name_of_particle_nr_get(tng
, ind
[i
], temp_name
, 256);
391 if (stat
!= TNG_SUCCESS
)
395 stat
= tng_atom_type_of_particle_nr_get(tng
, ind
[i
], temp_type
, 256);
396 if (stat
!= TNG_SUCCESS
)
400 tng_residue_atom_w_id_add(tng
, res
, temp_name
, temp_type
, ind
[i
], &atom
);
402 tng_molecule_existing_add(tng
, &mol
);
404 /* Set the count of the molecule containing the selected atoms to 1 and all
405 * other molecules to 0 */
406 tng_molecule_cnt_set(tng
, mol
, 1);
407 tng_num_molecule_types_get(tng
, &nMols
);
408 for (gmx_int64_t k
= 0; k
< nMols
; k
++)
410 tng_molecule_of_index_get(tng
, k
, &iterMol
);
415 tng_molecule_cnt_set(tng
, iterMol
, 0);
418 GMX_UNUSED_VALUE(tng
);
419 GMX_UNUSED_VALUE(nind
);
420 GMX_UNUSED_VALUE(ind
);
421 GMX_UNUSED_VALUE(name
);
425 /* TODO: If/when TNG acquires the ability to copy data blocks without
426 * uncompressing them, then this implemenation should be reconsidered.
427 * Ideally, gmx trjconv -f a.tng -o b.tng -b 10 -e 20 would be fast
428 * and lose no information. */
429 gmx_bool
gmx_read_next_tng_frame(tng_trajectory_t input
,
431 gmx_int64_t
*requestedIds
,
436 tng_function_status stat
;
437 gmx_int64_t numberOfAtoms
= -1, frameNumber
= -1;
438 gmx_int64_t nBlocks
, blockId
, *blockIds
= NULL
, codecId
;
441 double frameTime
= -1.0;
442 int size
, blockDependency
;
444 const int defaultNumIds
= 5;
445 static gmx_int64_t fallbackRequestedIds
[defaultNumIds
] =
447 TNG_TRAJ_BOX_SHAPE
, TNG_TRAJ_POSITIONS
,
448 TNG_TRAJ_VELOCITIES
, TNG_TRAJ_FORCES
,
463 /* If no specific IDs were requested read all block types that can
464 * currently be interpreted */
465 if (!requestedIds
|| numRequestedIds
== 0)
467 numRequestedIds
= defaultNumIds
;
468 requestedIds
= fallbackRequestedIds
;
471 stat
= tng_num_particles_get(input
, &numberOfAtoms
);
472 if (stat
!= TNG_SUCCESS
)
474 gmx_file("Cannot determine number of atoms from TNG file.");
476 fr
->natoms
= numberOfAtoms
;
478 if (!gmx_get_tng_data_block_types_of_next_frame(input
,
494 for (gmx_int64_t i
= 0; i
< nBlocks
; i
++)
496 blockId
= blockIds
[i
];
497 tng_data_block_dependency_get(input
, blockId
, &blockDependency
);
498 if (blockDependency
& TNG_PARTICLE_DEPENDENT
)
500 stat
= tng_util_particle_data_next_frame_read(input
,
509 stat
= tng_util_non_particle_data_next_frame_read(input
,
516 if (stat
== TNG_CRITICAL
)
518 gmx_file("Cannot read positions from TNG file.");
521 else if (stat
== TNG_FAILURE
)
527 case TNG_TRAJ_BOX_SHAPE
:
531 size
= sizeof(gmx_int64_t
);
534 size
= sizeof(float);
536 case TNG_DOUBLE_DATA
:
537 size
= sizeof(double);
540 gmx_incons("Illegal datatype of box shape values!");
542 for (int i
= 0; i
< DIM
; i
++)
544 convert_array_to_real_array(reinterpret_cast<char *>(values
) + size
* i
* DIM
,
545 reinterpret_cast<real
*>(fr
->box
[i
]),
546 getDistanceScaleFactor(input
),
553 case TNG_TRAJ_POSITIONS
:
554 srenew(fr
->x
, fr
->natoms
);
555 convert_array_to_real_array(values
,
556 reinterpret_cast<real
*>(fr
->x
),
557 getDistanceScaleFactor(input
),
562 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &prec
);
563 /* This must be updated if/when more lossy compression methods are added */
564 if (codecId
== TNG_TNG_COMPRESSION
)
570 case TNG_TRAJ_VELOCITIES
:
571 srenew(fr
->v
, fr
->natoms
);
572 convert_array_to_real_array(values
,
574 getDistanceScaleFactor(input
),
579 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &prec
);
580 /* This must be updated if/when more lossy compression methods are added */
581 if (codecId
== TNG_TNG_COMPRESSION
)
587 case TNG_TRAJ_FORCES
:
588 srenew(fr
->f
, fr
->natoms
);
589 convert_array_to_real_array(values
,
590 reinterpret_cast<real
*>(fr
->f
),
591 getDistanceScaleFactor(input
),
601 fr
->lambda
= *(reinterpret_cast<float *>(values
));
603 case TNG_DOUBLE_DATA
:
604 fr
->lambda
= *(reinterpret_cast<double *>(values
));
607 gmx_incons("Illegal datatype lambda value!");
612 gmx_warning("Illegal block type! Currently GROMACS tools can only handle certain data types. Skipping block.");
614 /* values does not have to be freed before reading next frame. It will
615 * be reallocated if it is not NULL. */
618 fr
->step
= static_cast<int>(frameNumber
);
620 // Convert the time to ps
621 fr
->time
= frameTime
/ PICO
;
624 /* values must be freed before leaving this function */
629 GMX_UNUSED_VALUE(input
);
630 GMX_UNUSED_VALUE(fr
);
631 GMX_UNUSED_VALUE(requestedIds
);
632 GMX_UNUSED_VALUE(numRequestedIds
);
637 void gmx_print_tng_molecule_system(tng_trajectory_t input
,
641 gmx_int64_t nMolecules
, nChains
, nResidues
, nAtoms
, *molCntList
;
642 tng_molecule_t molecule
;
644 tng_residue_t residue
;
646 char str
[256], varNAtoms
;
648 tng_num_molecule_types_get(input
, &nMolecules
);
649 tng_molecule_cnt_list_get(input
, &molCntList
);
650 /* Can the number of particles change in the trajectory or is it constant? */
651 tng_num_particles_variable_get(input
, &varNAtoms
);
653 for (gmx_int64_t i
= 0; i
< nMolecules
; i
++)
655 tng_molecule_of_index_get(input
, i
, &molecule
);
656 tng_molecule_name_get(input
, molecule
, str
, 256);
657 if (varNAtoms
== TNG_CONSTANT_N_ATOMS
)
659 if ((int)molCntList
[i
] == 0)
663 fprintf(stream
, "Molecule: %s, count: %d\n", str
, (int)molCntList
[i
]);
667 fprintf(stream
, "Molecule: %s\n", str
);
669 tng_molecule_num_chains_get(input
, molecule
, &nChains
);
672 for (gmx_int64_t j
= 0; j
< nChains
; j
++)
674 tng_molecule_chain_of_index_get(input
, molecule
, j
, &chain
);
675 tng_chain_name_get(input
, chain
, str
, 256);
676 fprintf(stream
, "\tChain: %s\n", str
);
677 tng_chain_num_residues_get(input
, chain
, &nResidues
);
678 for (gmx_int64_t k
= 0; k
< nResidues
; k
++)
680 tng_chain_residue_of_index_get(input
, chain
, k
, &residue
);
681 tng_residue_name_get(input
, residue
, str
, 256);
682 fprintf(stream
, "\t\tResidue: %s\n", str
);
683 tng_residue_num_atoms_get(input
, residue
, &nAtoms
);
684 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
686 tng_residue_atom_of_index_get(input
, residue
, l
, &atom
);
687 tng_atom_name_get(input
, atom
, str
, 256);
688 fprintf(stream
, "\t\t\tAtom: %s", str
);
689 tng_atom_type_get(input
, atom
, str
, 256);
690 fprintf(stream
, " (%s)\n", str
);
695 /* It is possible to have a molecule without chains, in which case
696 * residues in the molecule can be iterated through without going
700 tng_molecule_num_residues_get(input
, molecule
, &nResidues
);
703 for (gmx_int64_t k
= 0; k
< nResidues
; k
++)
705 tng_molecule_residue_of_index_get(input
, molecule
, k
, &residue
);
706 tng_residue_name_get(input
, residue
, str
, 256);
707 fprintf(stream
, "\t\tResidue: %s\n", str
);
708 tng_residue_num_atoms_get(input
, residue
, &nAtoms
);
709 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
711 tng_residue_atom_of_index_get(input
, residue
, l
, &atom
);
712 tng_atom_name_get(input
, atom
, str
, 256);
713 fprintf(stream
, "\t\t\tAtom: %s", str
);
714 tng_atom_type_get(input
, atom
, str
, 256);
715 fprintf(stream
, " (%s)\n", str
);
721 tng_molecule_num_atoms_get(input
, molecule
, &nAtoms
);
722 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
724 tng_molecule_atom_of_index_get(input
, molecule
, l
, &atom
);
725 tng_atom_name_get(input
, atom
, str
, 256);
726 fprintf(stream
, "\t\t\tAtom: %s", str
);
727 tng_atom_type_get(input
, atom
, str
, 256);
728 fprintf(stream
, " (%s)\n", str
);
734 GMX_UNUSED_VALUE(input
);
735 GMX_UNUSED_VALUE(stream
);
739 gmx_bool
gmx_get_tng_data_block_types_of_next_frame(tng_trajectory_t input
,
742 gmx_int64_t
*requestedIds
,
743 gmx_int64_t
*nextFrame
,
744 gmx_int64_t
*nBlocks
,
745 gmx_int64_t
**blockIds
)
748 tng_function_status stat
;
750 stat
= tng_util_trajectory_next_frame_present_data_blocks_find(input
, frame
,
751 nRequestedIds
, requestedIds
,
755 if (stat
== TNG_CRITICAL
)
757 gmx_file("Cannot read TNG file. Cannot find data blocks of next frame.");
759 else if (stat
== TNG_FAILURE
)
765 GMX_UNUSED_VALUE(input
);
766 GMX_UNUSED_VALUE(frame
);
767 GMX_UNUSED_VALUE(nRequestedIds
);
768 GMX_UNUSED_VALUE(requestedIds
);
769 GMX_UNUSED_VALUE(nextFrame
);
770 GMX_UNUSED_VALUE(nBlocks
);
771 GMX_UNUSED_VALUE(blockIds
);
776 gmx_bool
gmx_get_tng_data_next_frame_of_block_type(tng_trajectory_t input
,
779 gmx_int64_t
*frameNumber
,
781 gmx_int64_t
*nValuesPerFrame
,
789 tng_function_status stat
;
796 stat
= tng_data_block_name_get(input
, blockId
, name
, maxLen
);
797 if (stat
!= TNG_SUCCESS
)
799 gmx_file("Cannot read next frame of TNG file");
801 stat
= tng_data_block_dependency_get(input
, blockId
, &blockDependency
);
802 if (stat
!= TNG_SUCCESS
)
804 gmx_file("Cannot read next frame of TNG file");
806 if (blockDependency
& TNG_PARTICLE_DEPENDENT
)
808 tng_num_particles_get(input
, nAtoms
);
809 stat
= tng_util_particle_data_next_frame_read(input
,
818 *nAtoms
= 1; /* There are not actually any atoms, but it is used for
820 stat
= tng_util_non_particle_data_next_frame_read(input
,
827 if (stat
== TNG_CRITICAL
)
829 gmx_file("Cannot read next frame of TNG file");
831 if (stat
== TNG_FAILURE
)
837 stat
= tng_data_block_num_values_per_frame_get(input
, blockId
, nValuesPerFrame
);
838 if (stat
!= TNG_SUCCESS
)
840 gmx_file("Cannot read next frame of TNG file");
842 snew(*values
, sizeof(real
) * *nValuesPerFrame
* *nAtoms
);
843 convert_array_to_real_array(data
,
845 getDistanceScaleFactor(input
),
850 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &localPrec
);
852 /* This must be updated if/when more lossy compression methods are added */
853 if (codecId
!= TNG_TNG_COMPRESSION
)
865 GMX_UNUSED_VALUE(input
);
866 GMX_UNUSED_VALUE(blockId
);
867 GMX_UNUSED_VALUE(values
);
868 GMX_UNUSED_VALUE(frameNumber
);
869 GMX_UNUSED_VALUE(frameTime
);
870 GMX_UNUSED_VALUE(nValuesPerFrame
);
871 GMX_UNUSED_VALUE(nAtoms
);
872 GMX_UNUSED_VALUE(prec
);
873 GMX_UNUSED_VALUE(name
);
874 GMX_UNUSED_VALUE(maxLen
);
875 GMX_UNUSED_VALUE(bOK
);