2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2013,2014, 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.
35 #include "tngio_for_tools.h"
47 #include "../../external/tng_io/include/tng_io.h"
50 #include "gromacs/utility/common.h"
51 #include "gromacs/legacyheaders/types/atoms.h"
52 #include "gromacs/legacyheaders/smalloc.h"
53 #include "gromacs/legacyheaders/physics.h"
54 #include "gromacs/legacyheaders/gmx_fatal.h"
56 void gmx_prepare_tng_writing(const char *filename
,
58 tng_trajectory_t
*input
,
59 tng_trajectory_t
*output
,
61 const gmx_mtop_t
*mtop
,
63 const char *indexGroupName
)
66 /* FIXME after 5.0: Currently only standard block types are read */
67 const int defaultNumIds
= 5;
68 static gmx_int64_t fallbackIds
[defaultNumIds
] =
70 TNG_TRAJ_BOX_SHAPE
, TNG_TRAJ_POSITIONS
,
71 TNG_TRAJ_VELOCITIES
, TNG_TRAJ_FORCES
,
74 static char fallbackNames
[defaultNumIds
][32] =
76 "BOX SHAPE", "POSITIONS", "VELOCITIES",
81 gmx_tng_open(filename
, mode
, output
);
83 /* Do we have an input file in TNG format? If so, then there's
84 more data we can copy over, rather than having to improvise. */
87 /* Set parameters (compression, time per frame, molecule
88 * information, number of frames per frame set and writing
89 * intervals of positions, box shape and lambdas) of the
90 * output tng container based on their respective values int
91 * the input tng container */
92 double time
, compression_precision
;
93 gmx_int64_t n_frames_per_frame_set
, interval
= -1;
95 tng_compression_precision_get(*input
, &compression_precision
);
96 tng_compression_precision_set(*output
, compression_precision
);
97 // TODO make this configurable in a future version
98 char compression_type
= TNG_TNG_COMPRESSION
;
100 tng_molecule_system_copy(*input
, *output
);
102 tng_time_per_frame_get(*input
, &time
);
103 tng_time_per_frame_set(*output
, time
);
105 tng_num_frames_per_frame_set_get(*input
, &n_frames_per_frame_set
);
106 tng_num_frames_per_frame_set_set(*output
, n_frames_per_frame_set
);
108 for (int i
= 0; i
< defaultNumIds
; i
++)
110 if (tng_data_get_stride_length(*input
, fallbackIds
[i
], -1, &interval
)
113 switch (fallbackIds
[i
])
115 case TNG_TRAJ_POSITIONS
:
116 case TNG_TRAJ_VELOCITIES
:
117 tng_util_generic_write_interval_set(*output
, interval
, 3, fallbackIds
[i
],
118 fallbackNames
[i
], TNG_PARTICLE_BLOCK_DATA
,
121 case TNG_TRAJ_FORCES
:
122 tng_util_generic_write_interval_set(*output
, interval
, 3, fallbackIds
[i
],
123 fallbackNames
[i
], TNG_PARTICLE_BLOCK_DATA
,
124 TNG_GZIP_COMPRESSION
);
126 case TNG_TRAJ_BOX_SHAPE
:
127 tng_util_generic_write_interval_set(*output
, interval
, 9, fallbackIds
[i
],
128 fallbackNames
[i
], TNG_NON_PARTICLE_BLOCK_DATA
,
129 TNG_GZIP_COMPRESSION
);
132 tng_util_generic_write_interval_set(*output
, interval
, 1, fallbackIds
[i
],
133 fallbackNames
[i
], TNG_NON_PARTICLE_BLOCK_DATA
,
134 TNG_GZIP_COMPRESSION
);
144 /* TODO after trjconv is modularized: fix this so the user can
145 change precision when they are doing an operation where
146 this makes sense, and not otherwise.
148 char compression = bUseLossyCompression ? TNG_TNG_COMPRESSION : TNG_GZIP_COMPRESSION;
149 gmx_tng_set_compression_precision(*output, ndec2prec(nDecimalsOfPrecision));
151 gmx_tng_add_mtop(*output
, mtop
);
152 tng_num_frames_per_frame_set_set(*output
, 1);
155 if (index
&& nAtoms
> 0)
157 gmx_tng_setup_atom_subgroup(*output
, nAtoms
, index
, indexGroupName
);
160 /* If for some reason there are more requested atoms than there are atoms in the
161 * molecular system create a number of implicit atoms (without atom data) to
162 * compensate for that. */
165 tng_implicit_num_particles_set(*output
, nAtoms
);
168 GMX_UNUSED_VALUE(filename
);
169 GMX_UNUSED_VALUE(mode
);
170 GMX_UNUSED_VALUE(input
);
171 GMX_UNUSED_VALUE(output
);
172 GMX_UNUSED_VALUE(nAtoms
);
176 void gmx_write_tng_from_trxframe(tng_trajectory_t output
,
183 double timePerFrame
= frame
->time
* PICO
/ frame
->step
;
184 tng_time_per_frame_set(output
, timePerFrame
);
188 natoms
= frame
->natoms
;
190 gmx_fwrite_tng(output
,
195 (const rvec
*) frame
->box
,
197 (const rvec
*) frame
->x
,
198 (const rvec
*) frame
->v
,
199 (const rvec
*) frame
->f
);
201 GMX_UNUSED_VALUE(output
);
202 GMX_UNUSED_VALUE(frame
);
208 convert_array_to_real_array(void *from
,
220 if (sizeof(real
) == sizeof(float))
224 memcpy(to
, from
, nValues
* sizeof(real
) * nAtoms
);
228 for (i
= 0; i
< nAtoms
; i
++)
230 for (j
= 0; j
< nValues
; j
++)
232 to
[i
*nValues
+j
] = (real
)((float *)from
)[i
*nValues
+j
] * fact
;
239 for (i
= 0; i
< nAtoms
; i
++)
241 for (j
= 0; j
< nValues
; j
++)
243 to
[i
*nValues
+j
] = (real
)((float *)from
)[i
*nValues
+j
] * fact
;
249 for (i
= 0; i
< nAtoms
; i
++)
251 for (j
= 0; j
< nValues
; j
++)
253 to
[i
*nValues
+j
] = (real
)((gmx_int64_t
*)from
)[i
*nValues
+j
] * fact
;
257 case TNG_DOUBLE_DATA
:
258 if (sizeof(real
) == sizeof(double))
262 memcpy(to
, from
, nValues
* sizeof(real
) * nAtoms
);
266 for (i
= 0; i
< nAtoms
; i
++)
268 for (j
= 0; j
< nValues
; j
++)
270 to
[i
*nValues
+j
] = (real
)((double *)from
)[i
*nValues
+j
] * fact
;
277 for (i
= 0; i
< nAtoms
; i
++)
279 for (j
= 0; j
< nValues
; j
++)
281 to
[i
*nValues
+j
] = (real
)((double *)from
)[i
*nValues
+j
] * fact
;
287 gmx_incons("Illegal datatype when converting values to a real array!");
293 static real
getDistanceScaleFactor(tng_trajectory_t in
)
295 gmx_int64_t exp
= -1;
296 real distanceScaleFactor
;
298 // TODO Hopefully, TNG 2.0 will do this kind of thing for us
299 tng_distance_unit_exponential_get(in
, &exp
);
301 // GROMACS expects distances in nm
305 distanceScaleFactor
= NANO
/NANO
;
308 distanceScaleFactor
= NANO
/ANGSTROM
;
311 distanceScaleFactor
= pow(10.0, exp
+ 9.0);
314 return distanceScaleFactor
;
318 void gmx_tng_setup_atom_subgroup(tng_trajectory_t tng
,
324 gmx_int64_t nAtoms
, cnt
, nMols
;
325 tng_molecule_t mol
, iterMol
;
329 tng_function_status stat
;
331 tng_num_particles_get(tng
, &nAtoms
);
338 stat
= tng_molecule_find(tng
, name
, -1, &mol
);
339 if (stat
== TNG_SUCCESS
)
341 tng_molecule_num_atoms_get(tng
, mol
, &nAtoms
);
342 tng_molecule_cnt_get(tng
, mol
, &cnt
);
352 if (stat
== TNG_FAILURE
)
354 /* The indexed atoms are added to one separate molecule. */
355 tng_molecule_alloc(tng
, &mol
);
356 tng_molecule_name_set(tng
, mol
, name
);
357 tng_molecule_chain_add(tng
, mol
, "", &chain
);
359 for (int i
= 0; i
< nind
; i
++)
361 char temp_name
[256], temp_type
[256];
363 /* Try to retrieve the residue name of the atom */
364 stat
= tng_residue_name_of_particle_nr_get(tng
, ind
[i
], temp_name
, 256);
365 if (stat
!= TNG_SUCCESS
)
369 /* Check if the molecule of the selection already contains this residue */
370 if (tng_chain_residue_find(tng
, chain
, temp_name
, -1, &res
)
373 tng_chain_residue_add(tng
, chain
, temp_name
, &res
);
375 /* Try to find the original name and type of the atom */
376 stat
= tng_atom_name_of_particle_nr_get(tng
, ind
[i
], temp_name
, 256);
377 if (stat
!= TNG_SUCCESS
)
381 stat
= tng_atom_type_of_particle_nr_get(tng
, ind
[i
], temp_type
, 256);
382 if (stat
!= TNG_SUCCESS
)
386 tng_residue_atom_w_id_add(tng
, res
, temp_name
, temp_type
, ind
[i
], &atom
);
388 tng_molecule_existing_add(tng
, &mol
);
390 /* Set the count of the molecule containing the selected atoms to 1 and all
391 * other molecules to 0 */
392 tng_molecule_cnt_set(tng
, mol
, 1);
393 tng_num_molecule_types_get(tng
, &nMols
);
394 for (gmx_int64_t k
= 0; k
< nMols
; k
++)
396 tng_molecule_of_index_get(tng
, k
, &iterMol
);
401 tng_molecule_cnt_set(tng
, iterMol
, 0);
404 GMX_UNUSED_VALUE(tng
);
405 GMX_UNUSED_VALUE(nind
);
406 GMX_UNUSED_VALUE(ind
);
407 GMX_UNUSED_VALUE(name
);
411 /* TODO: If/when TNG acquires the ability to copy data blocks without
412 * uncompressing them, then this implemenation should be reconsidered.
413 * Ideally, gmx trjconv -f a.tng -o b.tng -b 10 -e 20 would be fast
414 * and lose no information. */
415 gmx_bool
gmx_read_next_tng_frame(tng_trajectory_t input
,
417 gmx_int64_t
*requestedIds
,
422 tng_function_status stat
;
423 gmx_int64_t numberOfAtoms
= -1, frameNumber
= -1;
424 gmx_int64_t nBlocks
, blockId
, *blockIds
= NULL
, codecId
;
427 double frameTime
= -1.0;
428 int size
, blockDependency
;
430 const int defaultNumIds
= 5;
431 static gmx_int64_t fallbackRequestedIds
[defaultNumIds
] =
433 TNG_TRAJ_BOX_SHAPE
, TNG_TRAJ_POSITIONS
,
434 TNG_TRAJ_VELOCITIES
, TNG_TRAJ_FORCES
,
449 /* If no specific IDs were requested read all block types that can
450 * currently be interpreted */
451 if (!requestedIds
|| numRequestedIds
== 0)
453 numRequestedIds
= defaultNumIds
;
454 requestedIds
= fallbackRequestedIds
;
457 stat
= tng_num_particles_get(input
, &numberOfAtoms
);
458 if (stat
!= TNG_SUCCESS
)
460 gmx_file("Cannot determine number of atoms from TNG file.");
462 fr
->natoms
= numberOfAtoms
;
464 if (!gmx_get_tng_data_block_types_of_next_frame(input
,
480 for (gmx_int64_t i
= 0; i
< nBlocks
; i
++)
482 blockId
= blockIds
[i
];
483 tng_data_block_dependency_get(input
, blockId
, &blockDependency
);
484 if (blockDependency
& TNG_PARTICLE_DEPENDENT
)
486 stat
= tng_util_particle_data_next_frame_read(input
,
495 stat
= tng_util_non_particle_data_next_frame_read(input
,
502 if (stat
== TNG_CRITICAL
)
504 gmx_file("Cannot read positions from TNG file.");
507 else if (stat
== TNG_FAILURE
)
513 case TNG_TRAJ_BOX_SHAPE
:
517 size
= sizeof(gmx_int64_t
);
520 size
= sizeof(float);
522 case TNG_DOUBLE_DATA
:
523 size
= sizeof(double);
526 size
= 0; /* Just to make the compiler happy. */
527 gmx_incons("Illegal datatype of box shape values!");
529 for (int i
= 0; i
< DIM
; i
++)
531 convert_array_to_real_array((char *)(values
) + size
* i
* DIM
,
533 getDistanceScaleFactor(input
),
540 case TNG_TRAJ_POSITIONS
:
541 srenew(fr
->x
, fr
->natoms
);
542 convert_array_to_real_array(values
,
544 getDistanceScaleFactor(input
),
549 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &prec
);
550 /* This must be updated if/when more lossy compression methods are added */
551 if (codecId
== TNG_TNG_COMPRESSION
)
557 case TNG_TRAJ_VELOCITIES
:
558 srenew(fr
->v
, fr
->natoms
);
559 convert_array_to_real_array(values
,
561 getDistanceScaleFactor(input
),
566 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &prec
);
567 /* This must be updated if/when more lossy compression methods are added */
568 if (codecId
== TNG_TNG_COMPRESSION
)
574 case TNG_TRAJ_FORCES
:
575 srenew(fr
->f
, fr
->natoms
);
576 convert_array_to_real_array(values
,
578 getDistanceScaleFactor(input
),
588 fr
->lambda
= (*(float *)values
);
590 case TNG_DOUBLE_DATA
:
591 fr
->lambda
= (*(double *)values
);
594 gmx_incons("Illegal datatype lambda value!");
599 gmx_warning("Illegal block type! Currently GROMACS tools can only handle certain data types. Skipping block.");
601 /* values does not have to be freed before reading next frame. It will
602 * be reallocated if it is not NULL. */
605 fr
->step
= (int) frameNumber
;
607 // Convert the time to ps
608 fr
->time
= frameTime
/ PICO
;
611 /* values must be freed before leaving this function */
616 GMX_UNUSED_VALUE(input
);
617 GMX_UNUSED_VALUE(fr
);
618 GMX_UNUSED_VALUE(requestedIds
);
623 void gmx_print_tng_molecule_system(tng_trajectory_t input
,
627 gmx_int64_t nMolecules
, nChains
, nResidues
, nAtoms
, *molCntList
;
628 tng_molecule_t molecule
;
630 tng_residue_t residue
;
632 char str
[256], varNAtoms
;
634 tng_num_molecule_types_get(input
, &nMolecules
);
635 tng_molecule_cnt_list_get(input
, &molCntList
);
636 /* Can the number of particles change in the trajectory or is it constant? */
637 tng_num_particles_variable_get(input
, &varNAtoms
);
639 for (gmx_int64_t i
= 0; i
< nMolecules
; i
++)
641 tng_molecule_of_index_get(input
, i
, &molecule
);
642 tng_molecule_name_get(input
, molecule
, str
, 256);
643 if (varNAtoms
== TNG_CONSTANT_N_ATOMS
)
645 if ((int)molCntList
[i
] == 0)
649 fprintf(stream
, "Molecule: %s, count: %d\n", str
, (int)molCntList
[i
]);
653 fprintf(stream
, "Molecule: %s\n", str
);
655 tng_molecule_num_chains_get(input
, molecule
, &nChains
);
658 for (gmx_int64_t j
= 0; j
< nChains
; j
++)
660 tng_molecule_chain_of_index_get(input
, molecule
, j
, &chain
);
661 tng_chain_name_get(input
, chain
, str
, 256);
662 fprintf(stream
, "\tChain: %s\n", str
);
663 tng_chain_num_residues_get(input
, chain
, &nResidues
);
664 for (gmx_int64_t k
= 0; k
< nResidues
; k
++)
666 tng_chain_residue_of_index_get(input
, chain
, k
, &residue
);
667 tng_residue_name_get(input
, residue
, str
, 256);
668 fprintf(stream
, "\t\tResidue: %s\n", str
);
669 tng_residue_num_atoms_get(input
, residue
, &nAtoms
);
670 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
672 tng_residue_atom_of_index_get(input
, residue
, l
, &atom
);
673 tng_atom_name_get(input
, atom
, str
, 256);
674 fprintf(stream
, "\t\t\tAtom: %s", str
);
675 tng_atom_type_get(input
, atom
, str
, 256);
676 fprintf(stream
, " (%s)\n", str
);
681 /* It is possible to have a molecule without chains, in which case
682 * residues in the molecule can be iterated through without going
686 tng_molecule_num_residues_get(input
, molecule
, &nResidues
);
689 for (gmx_int64_t k
= 0; k
< nResidues
; k
++)
691 tng_molecule_residue_of_index_get(input
, molecule
, k
, &residue
);
692 tng_residue_name_get(input
, residue
, str
, 256);
693 fprintf(stream
, "\t\tResidue: %s\n", str
);
694 tng_residue_num_atoms_get(input
, residue
, &nAtoms
);
695 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
697 tng_residue_atom_of_index_get(input
, residue
, l
, &atom
);
698 tng_atom_name_get(input
, atom
, str
, 256);
699 fprintf(stream
, "\t\t\tAtom: %s", str
);
700 tng_atom_type_get(input
, atom
, str
, 256);
701 fprintf(stream
, " (%s)\n", str
);
707 tng_molecule_num_atoms_get(input
, molecule
, &nAtoms
);
708 for (gmx_int64_t l
= 0; l
< nAtoms
; l
++)
710 tng_molecule_atom_of_index_get(input
, molecule
, l
, &atom
);
711 tng_atom_name_get(input
, atom
, str
, 256);
712 fprintf(stream
, "\t\t\tAtom: %s", str
);
713 tng_atom_type_get(input
, atom
, str
, 256);
714 fprintf(stream
, " (%s)\n", str
);
720 GMX_UNUSED_VALUE(input
);
721 GMX_UNUSED_VALUE(stream
);
725 gmx_bool
gmx_get_tng_data_block_types_of_next_frame(tng_trajectory_t input
,
728 gmx_int64_t
*requestedIds
,
729 gmx_int64_t
*nextFrame
,
730 gmx_int64_t
*nBlocks
,
731 gmx_int64_t
**blockIds
)
734 tng_function_status stat
;
736 stat
= tng_util_trajectory_next_frame_present_data_blocks_find(input
, frame
,
737 nRequestedIds
, requestedIds
,
741 if (stat
== TNG_CRITICAL
)
743 gmx_file("Cannot read TNG file. Cannot find data blocks of next frame.");
745 else if (stat
== TNG_FAILURE
)
751 GMX_UNUSED_VALUE(input
);
752 GMX_UNUSED_VALUE(frame
);
753 GMX_UNUSED_VALUE(nRequestedIds
);
754 GMX_UNUSED_VALUE(requestedIds
);
755 GMX_UNUSED_VALUE(nextFrame
);
756 GMX_UNUSED_VALUE(nBlocks
);
757 GMX_UNUSED_VALUE(blockIds
);
762 gmx_bool
gmx_get_tng_data_next_frame_of_block_type(tng_trajectory_t input
,
765 gmx_int64_t
*frameNumber
,
767 gmx_int64_t
*nValuesPerFrame
,
775 tng_function_status stat
;
782 stat
= tng_data_block_name_get(input
, blockId
, name
, maxLen
);
783 if (stat
!= TNG_SUCCESS
)
785 gmx_file("Cannot read next frame of TNG file");
787 stat
= tng_data_block_dependency_get(input
, blockId
, &blockDependency
);
788 if (stat
!= TNG_SUCCESS
)
790 gmx_file("Cannot read next frame of TNG file");
792 if (blockDependency
& TNG_PARTICLE_DEPENDENT
)
794 tng_num_particles_get(input
, nAtoms
);
795 stat
= tng_util_particle_data_next_frame_read(input
,
804 *nAtoms
= 1; /* There are not actually any atoms, but it is used for
806 stat
= tng_util_non_particle_data_next_frame_read(input
,
813 if (stat
== TNG_CRITICAL
)
815 gmx_file("Cannot read next frame of TNG file");
817 if (stat
== TNG_FAILURE
)
823 stat
= tng_data_block_num_values_per_frame_get(input
, blockId
, nValuesPerFrame
);
824 if (stat
!= TNG_SUCCESS
)
826 gmx_file("Cannot read next frame of TNG file");
828 snew(*values
, sizeof(real
) * *nValuesPerFrame
* *nAtoms
);
829 convert_array_to_real_array(data
,
831 getDistanceScaleFactor(input
),
836 tng_util_frame_current_compression_get(input
, blockId
, &codecId
, &localPrec
);
838 /* This must be updated if/when more lossy compression methods are added */
839 if (codecId
!= TNG_TNG_COMPRESSION
)
851 GMX_UNUSED_VALUE(input
);
852 GMX_UNUSED_VALUE(blockId
);
853 GMX_UNUSED_VALUE(values
);
854 GMX_UNUSED_VALUE(frameNumber
);
855 GMX_UNUSED_VALUE(frameTime
);
856 GMX_UNUSED_VALUE(nValuesPerFrame
);
857 GMX_UNUSED_VALUE(nAtoms
);
858 GMX_UNUSED_VALUE(prec
);
859 GMX_UNUSED_VALUE(name
);
860 GMX_UNUSED_VALUE(maxLen
);
861 GMX_UNUSED_VALUE(bOK
);