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.
46 #include "../../external/tng_io/include/tng_io.h"
49 #include "gromacs/legacyheaders/copyrite.h"
50 #include "gromacs/legacyheaders/gmx_fatal.h"
51 #include "gromacs/legacyheaders/main.h"
52 #include "gromacs/legacyheaders/physics.h"
53 #include "gromacs/math/utilities.h"
54 #include "gromacs/utility/common.h"
55 #include "gromacs/utility/gmxassert.h"
56 #include "gromacs/utility/programcontext.h"
59 static const char *modeToVerb(char mode
)
73 gmx_fatal(FARGS
, "Invalid file opening mode %c", mode
);
78 void gmx_tng_open(const char *filename
,
80 tng_trajectory_t
*tng
)
83 /* First check whether we have to make a backup,
84 * only for writing, not for read or append.
89 /* only make backups for normal gromacs */
90 make_backup(filename
);
94 /* tng must not be pointing at already allocated memory.
95 * Memory will be allocated by tng_util_trajectory_open() and must
96 * later on be freed by tng_util_trajectory_close(). */
97 if (TNG_SUCCESS
!= tng_util_trajectory_open(filename
, mode
, tng
))
99 /* TNG does return more than one degree of error, but there is
100 no use case for GROMACS handling the non-fatal errors
103 "%s while opening %s for %s",
104 gmx_strerror("file"),
109 if (mode
== 'w' || mode
== 'a')
111 /* FIXME in TNG: When adding data to the header, subsequent blocks might get
112 * overwritten. This could be solved by moving the first trajectory
113 * frame set(s) to the end of the file. Could that cause other problems,
114 * e.g. when continuing a simulation? */
116 gmx_gethostname(hostname
, 256);
119 tng_first_computer_name_set(*tng
, hostname
);
121 /* TODO: This should be implemented when the above fixme is done (adding data to
125 // tng_last_computer_name_set(*tng, hostname);
128 char programInfo
[256];
129 const char *precisionString
= "";
131 precisionString
= " (double precision)";
133 sprintf(programInfo
, "%.100s, %.128s%.24s",
134 gmx::getProgramContext().displayName(),
135 GromacsVersion(), precisionString
);
138 tng_first_program_name_set(*tng
, programInfo
);
140 /* TODO: This should be implemented when the above fixme is done (adding data to
144 // tng_last_program_name_set(*tng, programInfo);
149 getlogin_r(username
, 256);
152 tng_first_user_name_set(*tng
, username
);
154 /* TODO: This should be implemented when the above fixme is done (adding data to
158 // tng_last_user_name_set(*tng, username);
163 gmx_file("GROMACS was compiled without TNG support, cannot handle this file type");
164 GMX_UNUSED_VALUE(filename
);
165 GMX_UNUSED_VALUE(mode
);
166 GMX_UNUSED_VALUE(tng
);
170 void gmx_tng_close(tng_trajectory_t
*tng
)
172 /* We have to check that tng is set because
173 * tng_util_trajectory_close wants to return a NULL in it, and
174 * gives a fatal error if it is NULL. */
178 tng_util_trajectory_close(tng
);
181 GMX_UNUSED_VALUE(tng
);
186 static void addTngMoleculeFromTopology(tng_trajectory_t tng
,
187 const char *moleculeName
,
188 const t_atoms
*atoms
,
189 gmx_int64_t numMolecules
,
190 tng_molecule_t
*tngMol
)
192 if (tng_molecule_add(tng
, moleculeName
, tngMol
) != TNG_SUCCESS
)
194 gmx_file("Cannot add molecule to TNG molecular system.");
197 /* FIXME: The TNG atoms should contain mass and atomB info (for free
198 * energy calculations), i.e. in when it's available in TNG (2.0). */
199 for (int atomIt
= 0; atomIt
< atoms
->nr
; atomIt
++)
201 const t_atom
*at
= &atoms
->atom
[atomIt
];
202 /* FIXME: Currently the TNG API can only add atoms belonging to a
203 * residue and chain. Wait for TNG 2.0*/
206 const t_resinfo
*resInfo
= &atoms
->resinfo
[at
->resind
];
207 char chainName
[2] = {resInfo
->chainid
, 0};
208 tng_chain_t tngChain
= NULL
;
209 tng_residue_t tngRes
= NULL
;
210 tng_atom_t tngAtom
= NULL
;
212 if (tng_molecule_chain_find (tng
, *tngMol
, chainName
,
213 (gmx_int64_t
)-1, &tngChain
) !=
216 tng_molecule_chain_add (tng
, *tngMol
, chainName
,
220 /* FIXME: When TNG supports both residue index and residue
221 * number the latter should be used. Wait for TNG 2.0*/
222 if (tng_chain_residue_find(tng
, tngChain
, *resInfo
->name
,
223 at
->resind
+ 1, &tngRes
)
226 tng_chain_residue_add(tng
, tngChain
, *resInfo
->name
, &tngRes
);
228 tng_residue_atom_add(tng
, tngRes
, *(atoms
->atomname
[atomIt
]), *(atoms
->atomtype
[atomIt
]), &tngAtom
);
231 tng_molecule_cnt_set(tng
, *tngMol
, numMolecules
);
234 void gmx_tng_add_mtop(tng_trajectory_t tng
,
235 const gmx_mtop_t
*mtop
)
238 const t_ilist
*ilist
;
243 /* No topology information available to add. */
247 for (int molIt
= 0; molIt
< mtop
->nmolblock
; molIt
++)
249 tng_molecule_t tngMol
= NULL
;
250 const gmx_moltype_t
*molType
=
251 &mtop
->moltype
[mtop
->molblock
[molIt
].type
];
253 /* Add a molecule to the TNG trajectory with the same name as the
254 * current molecule. */
255 addTngMoleculeFromTopology(tng
,
258 mtop
->molblock
[molIt
].nmol
,
261 /* Bonds have to be deduced from interactions (constraints etc). Different
262 * interactions have different sets of parameters. */
263 /* Constraints are specified using two atoms */
264 for (i
= 0; i
< F_NRE
; i
++)
268 ilist
= &molType
->ilist
[i
];
272 while (j
< ilist
->nr
)
274 tng_molecule_bond_add(tng
, tngMol
, ilist
->iatoms
[j
], ilist
->iatoms
[j
+1], &tngBond
);
280 /* Settle is described using three atoms */
281 ilist
= &molType
->ilist
[F_SETTLE
];
285 while (j
< ilist
->nr
)
287 tng_molecule_bond_add(tng
, tngMol
, ilist
->iatoms
[j
], ilist
->iatoms
[j
+1], &tngBond
);
288 tng_molecule_bond_add(tng
, tngMol
, ilist
->iatoms
[j
], ilist
->iatoms
[j
+2], &tngBond
);
295 /*! \libinternal \brief Compute greatest common divisor of n1 and n2
296 * if they are positive.
298 * If only one of n1 and n2 is positive, then return it.
299 * If neither n1 or n2 is positive, then return -1. */
301 greatest_common_divisor_if_positive(int n1
, int n2
)
305 return (0 >= n2
) ? -1 : n2
;
312 /* We have a non-trivial greatest common divisor to compute. */
313 return gmx_greatest_common_divisor(n1
, n2
);
316 /* By default try to write 100 frames (of actual output) in each frame set.
317 * This number is the number of outputs of the most frequently written data
318 * type per frame set.
319 * TODO for 5.1: Verify that 100 frames per frame set is efficient for most
320 * setups regarding compression efficiency and compression time. Make this
321 * a hidden command-line option? */
322 const int defaultFramesPerFrameSet
= 100;
324 /*! \libinternal \brief Set the number of frames per frame
325 * set according to output intervals.
326 * The default is that 100 frames are written of the data
327 * that is written most often. */
328 static void tng_set_frames_per_frame_set(tng_trajectory_t tng
,
329 const gmx_bool bUseLossyCompression
,
330 const t_inputrec
*ir
)
334 /* Set the number of frames per frame set to contain at least
335 * defaultFramesPerFrameSet of the lowest common denominator of
336 * the writing interval of positions and velocities. */
337 /* FIXME after 5.0: consider nstenergy also? */
338 if (bUseLossyCompression
)
340 gcd
= ir
->nstxout_compressed
;
344 gcd
= greatest_common_divisor_if_positive(ir
->nstxout
, ir
->nstvout
);
345 gcd
= greatest_common_divisor_if_positive(gcd
, ir
->nstfout
);
352 tng_num_frames_per_frame_set_set(tng
, gcd
* defaultFramesPerFrameSet
);
355 /*! \libinternal \brief Set the data-writing intervals, and number of
356 * frames per frame set */
357 static void set_writing_intervals(tng_trajectory_t tng
,
358 const gmx_bool bUseLossyCompression
,
359 const t_inputrec
*ir
)
361 /* Define pointers to specific writing functions depending on if we
362 * write float or double data */
363 typedef tng_function_status (*set_writing_interval_func_pointer
)(tng_trajectory_t
,
371 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_double_set
;
373 set_writing_interval_func_pointer set_writing_interval
= tng_util_generic_write_interval_set
;
375 int xout
, vout
, fout
;
376 // int gcd = -1, lowest = -1;
379 tng_set_frames_per_frame_set(tng
, bUseLossyCompression
, ir
);
381 if (bUseLossyCompression
)
383 xout
= ir
->nstxout_compressed
;
386 compression
= TNG_TNG_COMPRESSION
;
393 compression
= TNG_GZIP_COMPRESSION
;
397 set_writing_interval(tng
, xout
, 3, TNG_TRAJ_POSITIONS
,
398 "POSITIONS", TNG_PARTICLE_BLOCK_DATA
,
400 /* The design of TNG makes it awkward to try to write a box
401 * with multiple periodicities, which might be co-prime. Since
402 * the use cases for the box with a frame consisting only of
403 * velocities seem low, for now we associate box writing with
404 * position writing. */
405 set_writing_interval(tng
, xout
, 9, TNG_TRAJ_BOX_SHAPE
,
406 "BOX SHAPE", TNG_NON_PARTICLE_BLOCK_DATA
,
407 TNG_GZIP_COMPRESSION
);
408 /* TODO: if/when we write energies to TNG also, reconsider how
409 * and when box information is written, because GROMACS
410 * behaviour pre-5.0 was to write the box with every
411 * trajectory frame and every energy frame, and probably
412 * people depend on this. */
414 /* TODO: If we need to write lambda values at steps when
415 * positions (or other data) are not also being written, then
416 * code in mdoutf.c will need to match however that is
418 set_writing_interval(tng
, xout
, 1, TNG_GMX_LAMBDA
,
419 "LAMBDAS", TNG_NON_PARTICLE_BLOCK_DATA
,
420 TNG_GZIP_COMPRESSION
);
422 /* FIXME: gcd and lowest currently not used. */
423 // gcd = greatest_common_divisor_if_positive(gcd, xout);
424 // if (lowest < 0 || xout < lowest)
431 set_writing_interval(tng
, ir
->nstvout
, 3, TNG_TRAJ_VELOCITIES
,
432 "VELOCITIES", TNG_PARTICLE_BLOCK_DATA
,
435 /* FIXME: gcd and lowest currently not used. */
436 // gcd = greatest_common_divisor_if_positive(gcd, vout);
437 // if (lowest < 0 || vout < lowest)
444 set_writing_interval(tng
, ir
->nstfout
, 3, TNG_TRAJ_FORCES
,
445 "FORCES", TNG_PARTICLE_BLOCK_DATA
,
446 TNG_GZIP_COMPRESSION
);
448 /* FIXME: gcd and lowest currently not used. */
449 // gcd = greatest_common_divisor_if_positive(gcd, fout);
450 // if (lowest < 0 || fout < lowest)
455 /* FIXME: See above. gcd interval for lambdas is disabled. */
458 // /* Lambdas written at an interval of the lowest common denominator
459 // * of other output */
460 // set_writing_interval(tng, gcd, 1, TNG_GMX_LAMBDA,
461 // "LAMBDAS", TNG_NON_PARTICLE_BLOCK_DATA,
462 // TNG_GZIP_COMPRESSION);
464 // if (gcd < lowest / 10)
466 // gmx_warning("The lowest common denominator of trajectory output is "
467 // "every %d step(s), whereas the shortest output interval "
468 // "is every %d steps.", gcd, lowest);
474 void gmx_tng_prepare_md_writing(tng_trajectory_t tng
,
475 const gmx_mtop_t
*mtop
,
476 const t_inputrec
*ir
)
479 gmx_tng_add_mtop(tng
, mtop
);
480 set_writing_intervals(tng
, FALSE
, ir
);
481 tng_time_per_frame_set(tng
, ir
->delta_t
* PICO
);
483 GMX_UNUSED_VALUE(tng
);
484 GMX_UNUSED_VALUE(mtop
);
485 GMX_UNUSED_VALUE(ir
);
490 /* Create a TNG molecule representing the selection groups
492 static void add_selection_groups(tng_trajectory_t tng
,
493 const gmx_mtop_t
*mtop
)
495 const gmx_moltype_t
*molType
;
496 const t_atoms
*atoms
;
498 const t_resinfo
*resInfo
;
499 const t_ilist
*ilist
;
500 int nAtoms
= 0, i
= 0, j
, molIt
, atomIt
, nameIndex
;
502 tng_molecule_t mol
, iterMol
;
510 /* The name of the TNG molecule containing the selection group is the
511 * same as the name of the selection group. */
512 nameIndex
= *mtop
->groups
.grps
[egcCompressedX
].nm_ind
;
513 groupName
= *mtop
->groups
.grpname
[nameIndex
];
515 tng_molecule_alloc(tng
, &mol
);
516 tng_molecule_name_set(tng
, mol
, groupName
);
517 tng_molecule_chain_add(tng
, mol
, "", &chain
);
518 for (molIt
= 0; molIt
< mtop
->nmoltype
; molIt
++)
520 molType
= &mtop
->moltype
[mtop
->molblock
[molIt
].type
];
522 atoms
= &molType
->atoms
;
524 for (j
= 0; j
< mtop
->molblock
[molIt
].nmol
; j
++)
526 bool bAtomsAdded
= FALSE
;
527 for (atomIt
= 0; atomIt
< atoms
->nr
; atomIt
++, i
++)
532 if (ggrpnr(&mtop
->groups
, egcCompressedX
, i
) != 0)
536 at
= &atoms
->atom
[atomIt
];
539 resInfo
= &atoms
->resinfo
[at
->resind
];
540 /* FIXME: When TNG supports both residue index and residue
541 * number the latter should be used. */
542 res_name
= *resInfo
->name
;
543 res_id
= at
->resind
+ 1;
547 res_name
= (char *)"";
550 if (tng_chain_residue_find(tng
, chain
, res_name
, res_id
, &res
)
553 /* Since there is ONE chain for selection groups do not keep the
554 * original residue IDs - otherwise there might be conflicts. */
555 tng_chain_residue_add(tng
, chain
, res_name
, &res
);
557 tng_residue_atom_w_id_add(tng
, res
, *(atoms
->atomname
[atomIt
]),
558 *(atoms
->atomtype
[atomIt
]),
559 atom_offset
+ atomIt
, &atom
);
566 for (int k
= 0; k
< F_NRE
; k
++)
570 ilist
= &molType
->ilist
[k
];
574 while (l
< ilist
->nr
)
577 atom1
= ilist
->iatoms
[l
] + atom_offset
;
578 atom2
= ilist
->iatoms
[l
+1] + atom_offset
;
579 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom1
) == 0 &&
580 ggrpnr(&mtop
->groups
, egcCompressedX
, atom2
) == 0)
582 tng_molecule_bond_add(tng
, mol
, ilist
->iatoms
[l
],
583 ilist
->iatoms
[l
+1], &tngBond
);
590 /* Settle is described using three atoms */
591 ilist
= &molType
->ilist
[F_SETTLE
];
595 while (l
< ilist
->nr
)
597 int atom1
, atom2
, atom3
;
598 atom1
= ilist
->iatoms
[l
] + atom_offset
;
599 atom2
= ilist
->iatoms
[l
+1] + atom_offset
;
600 atom3
= ilist
->iatoms
[l
+2] + atom_offset
;
601 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom1
) == 0)
603 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom2
) == 0)
605 tng_molecule_bond_add(tng
, mol
, atom1
,
608 if (ggrpnr(&mtop
->groups
, egcCompressedX
, atom3
) == 0)
610 tng_molecule_bond_add(tng
, mol
, atom1
,
618 atom_offset
+= atoms
->nr
;
623 tng_molecule_existing_add(tng
, &mol
);
624 tng_molecule_cnt_set(tng
, mol
, 1);
625 tng_num_molecule_types_get(tng
, &nMols
);
626 for (gmx_int64_t k
= 0; k
< nMols
; k
++)
628 tng_molecule_of_index_get(tng
, k
, &iterMol
);
633 tng_molecule_cnt_set(tng
, iterMol
, 0);
638 tng_molecule_free(tng
, &mol
);
643 void gmx_tng_set_compression_precision(tng_trajectory_t tng
,
647 tng_compression_precision_set(tng
, prec
);
649 GMX_UNUSED_VALUE(tng
);
650 GMX_UNUSED_VALUE(prec
);
654 void gmx_tng_prepare_low_prec_writing(tng_trajectory_t tng
,
655 const gmx_mtop_t
*mtop
,
656 const t_inputrec
*ir
)
659 gmx_tng_add_mtop(tng
, mtop
);
660 add_selection_groups(tng
, mtop
);
661 set_writing_intervals(tng
, TRUE
, ir
);
662 tng_time_per_frame_set(tng
, ir
->delta_t
* PICO
);
663 gmx_tng_set_compression_precision(tng
, ir
->x_compression_precision
);
665 GMX_UNUSED_VALUE(tng
);
666 GMX_UNUSED_VALUE(mtop
);
667 GMX_UNUSED_VALUE(ir
);
671 void gmx_fwrite_tng(tng_trajectory_t tng
,
672 const gmx_bool bUseLossyCompression
,
674 real elapsedPicoSeconds
,
683 typedef tng_function_status (*write_data_func_pointer
)(tng_trajectory_t
,
693 static write_data_func_pointer write_data
= tng_util_generic_with_time_double_write
;
695 static write_data_func_pointer write_data
= tng_util_generic_with_time_write
;
697 double elapsedSeconds
= elapsedPicoSeconds
* PICO
;
698 gmx_int64_t nParticles
;
704 /* This function might get called when the type of the
705 compressed trajectory is actually XTC. So we exit and move
710 tng_num_particles_get(tng
, &nParticles
);
711 if (nAtoms
!= (int)nParticles
)
713 tng_implicit_num_particles_set(tng
, nAtoms
);
716 if (bUseLossyCompression
)
718 compression
= TNG_TNG_COMPRESSION
;
722 compression
= TNG_GZIP_COMPRESSION
;
725 /* The writing is done using write_data, which writes float or double
726 * depending on the GROMACS compilation. */
729 GMX_ASSERT(box
, "Need a non-NULL box if positions are written");
731 if (write_data(tng
, step
, elapsedSeconds
,
732 reinterpret_cast<const real
*>(x
),
733 3, TNG_TRAJ_POSITIONS
, "POSITIONS",
734 TNG_PARTICLE_BLOCK_DATA
,
735 compression
) != TNG_SUCCESS
)
737 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
739 /* TNG-MF1 compression only compresses positions and velocities. Use lossless
740 * compression for box shape regardless of output mode */
741 if (write_data(tng
, step
, elapsedSeconds
,
742 reinterpret_cast<const real
*>(box
),
743 9, TNG_TRAJ_BOX_SHAPE
, "BOX SHAPE",
744 TNG_NON_PARTICLE_BLOCK_DATA
,
745 TNG_GZIP_COMPRESSION
) != TNG_SUCCESS
)
747 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
753 if (write_data(tng
, step
, elapsedSeconds
,
754 reinterpret_cast<const real
*>(v
),
755 3, TNG_TRAJ_VELOCITIES
, "VELOCITIES",
756 TNG_PARTICLE_BLOCK_DATA
,
757 compression
) != TNG_SUCCESS
)
759 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
765 /* TNG-MF1 compression only compresses positions and velocities. Use lossless
766 * compression for forces regardless of output mode */
767 if (write_data(tng
, step
, elapsedSeconds
,
768 reinterpret_cast<const real
*>(f
),
769 3, TNG_TRAJ_FORCES
, "FORCES",
770 TNG_PARTICLE_BLOCK_DATA
,
771 TNG_GZIP_COMPRESSION
) != TNG_SUCCESS
)
773 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
777 /* TNG-MF1 compression only compresses positions and velocities. Use lossless
778 * compression for lambdas regardless of output mode */
779 if (write_data(tng
, step
, elapsedSeconds
,
780 reinterpret_cast<const real
*>(&lambda
),
781 1, TNG_GMX_LAMBDA
, "LAMBDAS",
782 TNG_NON_PARTICLE_BLOCK_DATA
,
783 TNG_GZIP_COMPRESSION
) != TNG_SUCCESS
)
785 gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
788 GMX_UNUSED_VALUE(tng
);
789 GMX_UNUSED_VALUE(bUseLossyCompression
);
790 GMX_UNUSED_VALUE(step
);
791 GMX_UNUSED_VALUE(elapsedPicoSeconds
);
792 GMX_UNUSED_VALUE(lambda
);
793 GMX_UNUSED_VALUE(box
);
794 GMX_UNUSED_VALUE(nAtoms
);
801 void fflush_tng(tng_trajectory_t tng
)
808 tng_frame_set_premature_write(tng
, TNG_USE_HASH
);
810 GMX_UNUSED_VALUE(tng
);
814 float gmx_tng_get_time_of_final_frame(tng_trajectory_t tng
)
821 tng_num_frames_get(tng
, &nFrames
);
822 tng_util_time_of_frame_get(tng
, nFrames
- 1, &time
);
827 GMX_UNUSED_VALUE(tng
);