/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2013,2014, 2015 by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2018,2019,2020,2021, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include <algorithm>
#include <iterator>
#include <memory>
+#include <numeric>
#include <vector>
#if GMX_USE_TNG
-#include "tng/tng_io.h"
+# include "tng/tng_io.h"
#endif
#include "gromacs/math/units.h"
#include "gromacs/utility/sysinfo.h"
#include "gromacs/utility/unique_cptr.h"
+#if !GMX_USE_TNG
+using tng_trajectory_t = void*;
+#endif
+
/*! \brief Gromacs Wrapper around tng datatype
*
* This could in principle hold any GROMACS-specific requirements not yet
*/
struct gmx_tng_trajectory
{
- tng_trajectory_t tng; //!< Actual TNG handle (pointer)
- bool lastStepDataIsValid; //!< True if lastStep has been set
- std::int64_t lastStep; //!< Index/step used for last frame
- bool lastTimeDataIsValid; //!< True if lastTime has been set
- double lastTime; //!< Time of last frame (TNG unit is seconds)
- bool timePerFrameIsSet; //!< True if we have set the time per frame
+ tng_trajectory_t tng; //!< Actual TNG handle (pointer)
+ bool lastStepDataIsValid; //!< True if lastStep has been set
+ std::int64_t lastStep; //!< Index/step used for last frame
+ bool lastTimeDataIsValid; //!< True if lastTime has been set
+ double lastTime; //!< Time of last frame (TNG unit is seconds)
+ bool timePerFrameIsSet; //!< True if we have set the time per frame
+ int boxOutputInterval; //!< Number of steps between the output of box size
+ int lambdaOutputInterval; //!< Number of steps between the output of lambdas
};
-static const char *modeToVerb(char mode)
+#if GMX_USE_TNG
+static const char* modeToVerb(char mode)
{
- const char *p;
+ const char* p;
switch (mode)
{
- case 'r':
- p = "reading";
- break;
- case 'w':
- p = "writing";
- break;
- case 'a':
- p = "appending";
- break;
- default:
- gmx_fatal(FARGS, "Invalid file opening mode %c", mode);
- p = "";
- break;
+ case 'r': p = "reading"; break;
+ case 'w': p = "writing"; break;
+ case 'a': p = "appending"; break;
+ default: gmx_fatal(FARGS, "Invalid file opening mode %c", mode);
}
return p;
}
+#endif
-void gmx_tng_open(const char *filename,
- char mode,
- gmx_tng_trajectory_t *gmx_tng)
+void gmx_tng_open(const char* filename, char mode, gmx_tng_trajectory_t* gmx_tng)
{
#if GMX_USE_TNG
/* First check whether we have to make a backup,
(*gmx_tng)->lastStepDataIsValid = false;
(*gmx_tng)->lastTimeDataIsValid = false;
(*gmx_tng)->timePerFrameIsSet = false;
- tng_trajectory_t * tng = &(*gmx_tng)->tng;
+ tng_trajectory_t* tng = &(*gmx_tng)->tng;
/* tng must not be pointing at already allocated memory.
* Memory will be allocated by tng_util_trajectory_open() and must
/* TNG does return more than one degree of error, but there is
no use case for GROMACS handling the non-fatal errors
gracefully. */
- gmx_fatal(FARGS,
- "File I/O error while opening %s for %s",
- filename,
- modeToVerb(mode));
+ gmx_fatal(FARGS, "File I/O error while opening %s for %s", filename, modeToVerb(mode));
}
if (mode == 'w' || mode == 'a')
}
char programInfo[256];
- const char *precisionString = "";
-#if GMX_DOUBLE
+ const char* precisionString = "";
+# if GMX_DOUBLE
precisionString = " (double precision)";
-#endif
- sprintf(programInfo, "%.100s %.128s%.24s",
- gmx::getProgramContext().displayName(),
- gmx_version(), precisionString);
+# endif
+ sprintf(programInfo, "%.100s %.128s%.24s", gmx::getProgramContext().displayName(), gmx_version(), precisionString);
if (mode == 'w')
{
tng_first_program_name_set(*tng, programInfo);
#endif
}
-void gmx_tng_close(gmx_tng_trajectory_t *gmx_tng)
+void gmx_tng_close(gmx_tng_trajectory_t* gmx_tng)
{
/* We have to check that tng is set because
* tng_util_trajectory_close wants to return a NULL in it, and
{
return;
}
- tng_trajectory_t * tng = &(*gmx_tng)->tng;
+ tng_trajectory_t* tng = &(*gmx_tng)->tng;
if (tng)
{
#if GMX_USE_TNG
static void addTngMoleculeFromTopology(gmx_tng_trajectory_t gmx_tng,
- const char *moleculeName,
- const t_atoms *atoms,
- gmx_int64_t numMolecules,
- tng_molecule_t *tngMol)
+ const char* moleculeName,
+ const t_atoms* atoms,
+ int64_t numMolecules,
+ tng_molecule_t* tngMol)
{
tng_trajectory_t tng = gmx_tng->tng;
tng_chain_t tngChain = nullptr;
for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++)
{
- const t_atom *at = &atoms->atom[atomIndex];
+ const t_atom* at = &atoms->atom[atomIndex];
/* FIXME: Currently the TNG API can only add atoms belonging to a
* residue and chain. Wait for TNG 2.0*/
if (atoms->nres > 0)
{
- const t_resinfo *resInfo = &atoms->resinfo[at->resind];
- char chainName[2] = {resInfo->chainid, 0};
- tng_atom_t tngAtom = nullptr;
- t_atom *prevAtom;
+ const t_resinfo* resInfo = &atoms->resinfo[at->resind];
+ char chainName[2] = { resInfo->chainid, 0 };
+ tng_atom_t tngAtom = nullptr;
+ t_atom* prevAtom;
if (atomIndex > 0)
{
{
/* If this is the first atom or if the chain changed add
* the chain to the TNG molecular system. */
- if (!prevAtom || resInfo->chainid !=
- atoms->resinfo[prevAtom->resind].chainid)
+ if (!prevAtom || resInfo->chainid != atoms->resinfo[prevAtom->resind].chainid)
{
- tng_molecule_chain_add(tng, *tngMol, chainName,
- &tngChain);
+ tng_molecule_chain_add(tng, *tngMol, chainName, &tngChain);
}
/* FIXME: When TNG supports both residue index and residue
* number the latter should be used. Wait for TNG 2.0*/
tng_chain_residue_add(tng, tngChain, *resInfo->name, &tngRes);
}
- tng_residue_atom_add(tng, tngRes, *(atoms->atomname[atomIndex]), *(atoms->atomtype[atomIndex]), &tngAtom);
+ tng_residue_atom_add(
+ tng, tngRes, *(atoms->atomname[atomIndex]), *(atoms->atomtype[atomIndex]), &tngAtom);
}
}
tng_molecule_cnt_set(tng, *tngMol, numMolecules);
}
-void gmx_tng_add_mtop(gmx_tng_trajectory_t gmx_tng,
- const gmx_mtop_t *mtop)
+void gmx_tng_add_mtop(gmx_tng_trajectory_t gmx_tng, const gmx_mtop_t* mtop)
{
- int i;
- int j;
- std::vector<real> atomCharges;
- std::vector<real> atomMasses;
- const t_ilist *ilist;
- tng_bond_t tngBond;
- char datatype;
+ int i;
+ int j;
+ std::vector<real> atomCharges;
+ std::vector<real> atomMasses;
+ tng_bond_t tngBond;
+ char datatype;
- tng_trajectory_t tng = gmx_tng->tng;
+ tng_trajectory_t tng = gmx_tng->tng;
if (!mtop)
{
return;
}
-#if GMX_DOUBLE
+# if GMX_DOUBLE
datatype = TNG_DOUBLE_DATA;
-#else
- datatype = TNG_FLOAT_DATA;
-#endif
+# else
+ datatype = TNG_FLOAT_DATA;
+# endif
atomCharges.reserve(mtop->natoms);
atomMasses.reserve(mtop->natoms);
- for (int molIndex = 0; molIndex < mtop->nmolblock; molIndex++)
+ for (const gmx_molblock_t& molBlock : mtop->molblock)
{
tng_molecule_t tngMol = nullptr;
- const gmx_moltype_t *molType = &mtop->moltype[mtop->molblock[molIndex].type];
+ const gmx_moltype_t* molType = &mtop->moltype[molBlock.type];
/* Add a molecule to the TNG trajectory with the same name as the
* current molecule. */
- addTngMoleculeFromTopology(gmx_tng,
- *(molType->name),
- &molType->atoms,
- mtop->molblock[molIndex].nmol,
- &tngMol);
+ addTngMoleculeFromTopology(gmx_tng, *(molType->name), &molType->atoms, molBlock.nmol, &tngMol);
/* Bonds have to be deduced from interactions (constraints etc). Different
* interactions have different sets of parameters. */
{
if (IS_CHEMBOND(i))
{
- ilist = &molType->ilist[i];
- if (ilist)
+ const InteractionList& ilist = molType->ilist[i];
+ j = 1;
+ while (j < ilist.size())
{
- j = 1;
- while (j < ilist->nr)
- {
- tng_molecule_bond_add(tng, tngMol, ilist->iatoms[j], ilist->iatoms[j+1], &tngBond);
- j += 3;
- }
+ tng_molecule_bond_add(tng, tngMol, ilist.iatoms[j], ilist.iatoms[j + 1], &tngBond);
+ j += 3;
}
}
}
/* Settle is described using three atoms */
- ilist = &molType->ilist[F_SETTLE];
- if (ilist)
+ const InteractionList& ilist = molType->ilist[F_SETTLE];
+ j = 1;
+ while (j < ilist.size())
{
- j = 1;
- while (j < ilist->nr)
- {
- tng_molecule_bond_add(tng, tngMol, ilist->iatoms[j], ilist->iatoms[j+1], &tngBond);
- tng_molecule_bond_add(tng, tngMol, ilist->iatoms[j], ilist->iatoms[j+2], &tngBond);
- j += 4;
- }
+ tng_molecule_bond_add(tng, tngMol, ilist.iatoms[j], ilist.iatoms[j + 1], &tngBond);
+ tng_molecule_bond_add(tng, tngMol, ilist.iatoms[j], ilist.iatoms[j + 2], &tngBond);
+ j += 4;
}
/* First copy atom charges and masses, first atom by atom and then copy the memory for the molecule instances.
* FIXME: Atom B state data should also be written to TNG (v 2.0?) */
atomCharges.push_back(molType->atoms.atom[atomCounter].q);
atomMasses.push_back(molType->atoms.atom[atomCounter].m);
}
- for (int molCounter = 1; molCounter < mtop->molblock[molIndex].nmol; molCounter++)
+ for (int molCounter = 1; molCounter < molBlock.nmol; molCounter++)
{
- std::copy_n(atomCharges.end() - molType->atoms.nr, molType->atoms.nr, std::back_inserter(atomCharges));
- std::copy_n(atomMasses.end() - molType->atoms.nr, molType->atoms.nr, std::back_inserter(atomMasses));
+ std::copy_n(atomCharges.end() - molType->atoms.nr,
+ molType->atoms.nr,
+ std::back_inserter(atomCharges));
+ std::copy_n(atomMasses.end() - molType->atoms.nr, molType->atoms.nr, std::back_inserter(atomMasses));
}
}
/* Write the TNG data blocks. */
- tng_particle_data_block_add(tng, TNG_TRAJ_PARTIAL_CHARGES, "PARTIAL CHARGES",
- datatype, TNG_NON_TRAJECTORY_BLOCK,
- 1, 1, 1, 0, mtop->natoms,
- TNG_GZIP_COMPRESSION, atomCharges.data());
- tng_particle_data_block_add(tng, TNG_TRAJ_MASSES, "ATOM MASSES",
- datatype, TNG_NON_TRAJECTORY_BLOCK,
- 1, 1, 1, 0, mtop->natoms,
- TNG_GZIP_COMPRESSION, atomMasses.data());
+ tng_particle_data_block_add(tng,
+ TNG_TRAJ_PARTIAL_CHARGES,
+ "PARTIAL CHARGES",
+ datatype,
+ TNG_NON_TRAJECTORY_BLOCK,
+ 1,
+ 1,
+ 1,
+ 0,
+ mtop->natoms,
+ TNG_GZIP_COMPRESSION,
+ atomCharges.data());
+ tng_particle_data_block_add(tng,
+ TNG_TRAJ_MASSES,
+ "ATOM MASSES",
+ datatype,
+ TNG_NON_TRAJECTORY_BLOCK,
+ 1,
+ 1,
+ 1,
+ 0,
+ mtop->natoms,
+ TNG_GZIP_COMPRESSION,
+ atomMasses.data());
}
/*! \libinternal \brief Compute greatest common divisor of n1 and n2
*
* If only one of n1 and n2 is positive, then return it.
* If neither n1 or n2 is positive, then return -1. */
-static int
-greatest_common_divisor_if_positive(int n1, int n2)
+static int greatest_common_divisor_if_positive(int n1, int n2)
{
if (0 >= n1)
{
}
/* We have a non-trivial greatest common divisor to compute. */
- return gmx_greatest_common_divisor(n1, n2);
+ return std::gcd(n1, n2);
}
/* By default try to write 100 frames (of actual output) in each frame set.
* set according to output intervals.
* The default is that 100 frames are written of the data
* that is written most often. */
-static void tng_set_frames_per_frame_set(gmx_tng_trajectory_t gmx_tng,
- const gmx_bool bUseLossyCompression,
- const t_inputrec *ir)
+static void tng_set_frames_per_frame_set(gmx_tng_trajectory_t gmx_tng,
+ const gmx_bool bUseLossyCompression,
+ const t_inputrec* ir)
{
int gcd = -1;
tng_trajectory_t tng = gmx_tng->tng;
/*! \libinternal \brief Set the data-writing intervals, and number of
* frames per frame set */
-static void set_writing_intervals(gmx_tng_trajectory_t gmx_tng,
- const gmx_bool bUseLossyCompression,
- const t_inputrec *ir)
+static void set_writing_intervals(gmx_tng_trajectory_t gmx_tng,
+ const gmx_bool bUseLossyCompression,
+ const t_inputrec* ir)
{
tng_trajectory_t tng = gmx_tng->tng;
/* Define pointers to specific writing functions depending on if we
* write float or double data */
- typedef tng_function_status (*set_writing_interval_func_pointer)(tng_trajectory_t,
- const gmx_int64_t,
- const gmx_int64_t,
- const gmx_int64_t,
- const char*,
- const char,
- const char);
-#if GMX_DOUBLE
+ typedef tng_function_status (*set_writing_interval_func_pointer)(
+ tng_trajectory_t, const int64_t, const int64_t, const int64_t, const char*, const char, const char);
+# if GMX_DOUBLE
set_writing_interval_func_pointer set_writing_interval = tng_util_generic_write_interval_double_set;
-#else
+# else
set_writing_interval_func_pointer set_writing_interval = tng_util_generic_write_interval_set;
-#endif
+# endif
int xout, vout, fout;
int gcd = -1, lowest = -1;
char compression;
if (bUseLossyCompression)
{
- xout = ir->nstxout_compressed;
+ xout = ir->nstxout_compressed;
/* If there is no uncompressed coordinate output write forces
and velocities to the compressed tng file. */
if (ir->nstxout)
{
- vout = 0;
- fout = 0;
+ vout = 0;
+ fout = 0;
}
else
{
- vout = ir->nstvout;
- fout = ir->nstfout;
+ vout = ir->nstvout;
+ fout = ir->nstfout;
}
compression = TNG_TNG_COMPRESSION;
}
}
if (xout)
{
- set_writing_interval(tng, xout, 3, TNG_TRAJ_POSITIONS,
- "POSITIONS", TNG_PARTICLE_BLOCK_DATA,
- compression);
+ set_writing_interval(
+ tng, xout, 3, TNG_TRAJ_POSITIONS, "POSITIONS", TNG_PARTICLE_BLOCK_DATA, compression);
/* TODO: if/when we write energies to TNG also, reconsider how
* and when box information is written, because GROMACS
* behaviour pre-5.0 was to write the box with every
}
if (vout)
{
- set_writing_interval(tng, vout, 3, TNG_TRAJ_VELOCITIES,
- "VELOCITIES", TNG_PARTICLE_BLOCK_DATA,
- compression);
+ set_writing_interval(
+ tng, vout, 3, TNG_TRAJ_VELOCITIES, "VELOCITIES", TNG_PARTICLE_BLOCK_DATA, compression);
gcd = greatest_common_divisor_if_positive(gcd, vout);
if (lowest < 0 || vout < lowest)
}
if (fout)
{
- set_writing_interval(tng, fout, 3, TNG_TRAJ_FORCES,
- "FORCES", TNG_PARTICLE_BLOCK_DATA,
- TNG_GZIP_COMPRESSION);
+ set_writing_interval(
+ tng, fout, 3, TNG_TRAJ_FORCES, "FORCES", TNG_PARTICLE_BLOCK_DATA, TNG_GZIP_COMPRESSION);
gcd = greatest_common_divisor_if_positive(gcd, fout);
if (lowest < 0 || fout < lowest)
{
/* Lambdas and box shape written at an interval of the lowest common
denominator of other output */
- set_writing_interval(tng, gcd, 1, TNG_GMX_LAMBDA,
- "LAMBDAS", TNG_NON_PARTICLE_BLOCK_DATA,
- TNG_GZIP_COMPRESSION);
+ set_writing_interval(
+ tng, gcd, 1, TNG_GMX_LAMBDA, "LAMBDAS", TNG_NON_PARTICLE_BLOCK_DATA, TNG_GZIP_COMPRESSION);
- set_writing_interval(tng, gcd, 9, TNG_TRAJ_BOX_SHAPE,
- "BOX SHAPE", TNG_NON_PARTICLE_BLOCK_DATA,
- TNG_GZIP_COMPRESSION);
+ set_writing_interval(
+ tng, gcd, 9, TNG_TRAJ_BOX_SHAPE, "BOX SHAPE", TNG_NON_PARTICLE_BLOCK_DATA, TNG_GZIP_COMPRESSION);
+ gmx_tng->lambdaOutputInterval = gcd;
+ gmx_tng->boxOutputInterval = gcd;
if (gcd < lowest / 10)
{
- gmx_warning("The lowest common denominator of trajectory output is "
- "every %d step(s), whereas the shortest output interval "
- "is every %d steps.", gcd, lowest);
+ gmx_warning(
+ "The lowest common denominator of trajectory output is "
+ "every %d step(s), whereas the shortest output interval "
+ "is every %d steps.",
+ gcd,
+ lowest);
}
}
}
#endif
-void gmx_tng_prepare_md_writing(gmx_tng_trajectory_t gmx_tng,
- const gmx_mtop_t *mtop,
- const t_inputrec *ir)
+void gmx_tng_prepare_md_writing(gmx_tng_trajectory_t gmx_tng, const gmx_mtop_t* mtop, const t_inputrec* ir)
{
#if GMX_USE_TNG
gmx_tng_add_mtop(gmx_tng, mtop);
set_writing_intervals(gmx_tng, FALSE, ir);
- tng_time_per_frame_set(gmx_tng->tng, ir->delta_t * PICO);
+ tng_time_per_frame_set(gmx_tng->tng, ir->delta_t * gmx::c_pico);
gmx_tng->timePerFrameIsSet = true;
#else
GMX_UNUSED_VALUE(gmx_tng);
#if GMX_USE_TNG
/* Check if all atoms in the molecule system are selected
* by a selection group of type specified by gtype. */
-static gmx_bool all_atoms_selected(const gmx_mtop_t *mtop,
- const int gtype)
+static gmx_bool all_atoms_selected(const gmx_mtop_t* mtop, const SimulationAtomGroupType gtype)
{
- const gmx_moltype_t *molType;
- const t_atoms *atoms;
-
/* Iterate through all atoms in the molecule system and
* check if they belong to a selection group of the
* requested type. */
- for (int molIndex = 0, i = 0; molIndex < mtop->nmoltype; molIndex++)
+ int i = 0;
+ for (const gmx_molblock_t& molBlock : mtop->molblock)
{
- molType = &mtop->moltype[mtop->molblock[molIndex].type];
-
- atoms = &molType->atoms;
+ const gmx_moltype_t& molType = mtop->moltype[molBlock.type];
+ const t_atoms& atoms = molType.atoms;
- for (int j = 0; j < mtop->molblock[molIndex].nmol; j++)
+ for (int j = 0; j < molBlock.nmol; j++)
{
- for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++, i++)
+ for (int atomIndex = 0; atomIndex < atoms.nr; atomIndex++, i++)
{
- if (ggrpnr(&mtop->groups, gtype, i) != 0)
+ if (getGroupType(mtop->groups, gtype, i) != 0)
{
return FALSE;
}
* is egcCompressedX, but other selections should be added when
* e.g. writing energies is implemented.
*/
-static void add_selection_groups(gmx_tng_trajectory_t gmx_tng,
- const gmx_mtop_t *mtop)
+static void add_selection_groups(gmx_tng_trajectory_t gmx_tng, const gmx_mtop_t* mtop)
{
- const gmx_moltype_t *molType;
- const t_atoms *atoms;
- const t_atom *at;
- const t_resinfo *resInfo;
- const t_ilist *ilist;
- int nameIndex;
- int atom_offset = 0;
- tng_molecule_t mol, iterMol;
- tng_chain_t chain;
- tng_residue_t res;
- tng_atom_t atom;
- tng_bond_t tngBond;
- gmx_int64_t nMols;
- char *groupName;
- tng_trajectory_t tng = gmx_tng->tng;
+ const t_atoms* atoms;
+ const t_atom* at;
+ const t_resinfo* resInfo;
+ int nameIndex;
+ int atom_offset = 0;
+ tng_molecule_t mol, iterMol;
+ tng_chain_t chain;
+ tng_residue_t res;
+ tng_atom_t atom;
+ tng_bond_t tngBond;
+ int64_t nMols;
+ char* groupName;
+ tng_trajectory_t tng = gmx_tng->tng;
/* TODO: When the TNG molecules block is more flexible TNG selection
* groups should not need all atoms specified. It should be possible
/* If no atoms are selected we do not need to create a
* TNG selection group molecule. */
- if (mtop->groups.ngrpnr[egcCompressedX] == 0)
+ if (mtop->groups.numberOfGroupNumbers(SimulationAtomGroupType::CompressedPositionOutput) == 0)
{
return;
}
/* If all atoms are selected we do not have to create a selection
* group molecule in the TNG molecule system. */
- if (all_atoms_selected(mtop, egcCompressedX))
+ if (all_atoms_selected(mtop, SimulationAtomGroupType::CompressedPositionOutput))
{
return;
}
/* The name of the TNG molecule containing the selection group is the
* same as the name of the selection group. */
- nameIndex = *mtop->groups.grps[egcCompressedX].nm_ind;
- groupName = *mtop->groups.grpname[nameIndex];
+ nameIndex = mtop->groups.groups[SimulationAtomGroupType::CompressedPositionOutput][0];
+ groupName = *mtop->groups.groupNames[nameIndex];
tng_molecule_alloc(tng, &mol);
tng_molecule_name_set(tng, mol, groupName);
tng_molecule_chain_add(tng, mol, "", &chain);
- for (int molIndex = 0, i = 0; molIndex < mtop->nmoltype; molIndex++)
+ int i = 0;
+ for (const gmx_molblock_t& molBlock : mtop->molblock)
{
- molType = &mtop->moltype[mtop->molblock[molIndex].type];
+ const gmx_moltype_t& molType = mtop->moltype[molBlock.type];
- atoms = &molType->atoms;
+ atoms = &molType.atoms;
- for (int j = 0; j < mtop->molblock[molIndex].nmol; j++)
+ for (int j = 0; j < molBlock.nmol; j++)
{
bool bAtomsAdded = FALSE;
for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++, i++)
{
- char *res_name;
+ char* res_name;
int res_id;
- if (ggrpnr(&mtop->groups, egcCompressedX, i) != 0)
+ if (getGroupType(mtop->groups, SimulationAtomGroupType::CompressedPositionOutput, i) != 0)
{
continue;
}
}
else
{
- res_name = (char *)"";
+ res_name = const_cast<char*>("");
res_id = 0;
}
- if (tng_chain_residue_find(tng, chain, res_name, res_id, &res)
- != TNG_SUCCESS)
+ if (tng_chain_residue_find(tng, chain, res_name, res_id, &res) != TNG_SUCCESS)
{
/* Since there is ONE chain for selection groups do not keep the
* original residue IDs - otherwise there might be conflicts. */
tng_chain_residue_add(tng, chain, res_name, &res);
}
- tng_residue_atom_w_id_add(tng, res, *(atoms->atomname[atomIndex]),
+ tng_residue_atom_w_id_add(tng,
+ res,
+ *(atoms->atomname[atomIndex]),
*(atoms->atomtype[atomIndex]),
- atom_offset + atomIndex, &atom);
+ atom_offset + atomIndex,
+ &atom);
bAtomsAdded = TRUE;
}
/* Add bonds. */
{
if (IS_CHEMBOND(k))
{
- ilist = &molType->ilist[k];
- if (ilist)
+ const InteractionList& ilist = molType.ilist[k];
+ for (int l = 1; l < ilist.size(); l += 3)
{
- int l = 1;
- while (l < ilist->nr)
+ int atom1, atom2;
+ atom1 = ilist.iatoms[l] + atom_offset;
+ atom2 = ilist.iatoms[l + 1] + atom_offset;
+ if (getGroupType(mtop->groups, SimulationAtomGroupType::CompressedPositionOutput, atom1)
+ == 0
+ && getGroupType(mtop->groups, SimulationAtomGroupType::CompressedPositionOutput, atom2)
+ == 0)
{
- int atom1, atom2;
- atom1 = ilist->iatoms[l] + atom_offset;
- atom2 = ilist->iatoms[l+1] + atom_offset;
- if (ggrpnr(&mtop->groups, egcCompressedX, atom1) == 0 &&
- ggrpnr(&mtop->groups, egcCompressedX, atom2) == 0)
- {
- tng_molecule_bond_add(tng, mol, ilist->iatoms[l],
- ilist->iatoms[l+1], &tngBond);
- }
- l += 3;
+ tng_molecule_bond_add(
+ tng, mol, ilist.iatoms[l], ilist.iatoms[l + 1], &tngBond);
}
}
}
}
/* Settle is described using three atoms */
- ilist = &molType->ilist[F_SETTLE];
- if (ilist)
+ const InteractionList& ilist = molType.ilist[F_SETTLE];
+ for (int l = 1; l < ilist.size(); l += 4)
{
- int l = 1;
- while (l < ilist->nr)
+ int atom1, atom2, atom3;
+ atom1 = ilist.iatoms[l] + atom_offset;
+ atom2 = ilist.iatoms[l + 1] + atom_offset;
+ atom3 = ilist.iatoms[l + 2] + atom_offset;
+ if (getGroupType(mtop->groups, SimulationAtomGroupType::CompressedPositionOutput, atom1)
+ == 0)
{
- int atom1, atom2, atom3;
- atom1 = ilist->iatoms[l] + atom_offset;
- atom2 = ilist->iatoms[l+1] + atom_offset;
- atom3 = ilist->iatoms[l+2] + atom_offset;
- if (ggrpnr(&mtop->groups, egcCompressedX, atom1) == 0)
+ if (getGroupType(mtop->groups, SimulationAtomGroupType::CompressedPositionOutput, atom2)
+ == 0)
{
- if (ggrpnr(&mtop->groups, egcCompressedX, atom2) == 0)
- {
- tng_molecule_bond_add(tng, mol, atom1,
- atom2, &tngBond);
- }
- if (ggrpnr(&mtop->groups, egcCompressedX, atom3) == 0)
- {
- tng_molecule_bond_add(tng, mol, atom1,
- atom3, &tngBond);
- }
+ tng_molecule_bond_add(tng, mol, atom1, atom2, &tngBond);
+ }
+ if (getGroupType(mtop->groups, SimulationAtomGroupType::CompressedPositionOutput, atom3)
+ == 0)
+ {
+ tng_molecule_bond_add(tng, mol, atom1, atom3, &tngBond);
}
- l += 4;
}
}
}
tng_molecule_existing_add(tng, &mol);
tng_molecule_cnt_set(tng, mol, 1);
tng_num_molecule_types_get(tng, &nMols);
- for (gmx_int64_t k = 0; k < nMols; k++)
+ for (int64_t k = 0; k < nMols; k++)
{
tng_molecule_of_index_get(tng, k, &iterMol);
if (iterMol == mol)
}
#endif
-void gmx_tng_set_compression_precision(gmx_tng_trajectory_t gmx_tng,
- real prec)
+void gmx_tng_set_compression_precision(gmx_tng_trajectory_t gmx_tng, real prec)
{
#if GMX_USE_TNG
tng_compression_precision_set(gmx_tng->tng, prec);
#endif
}
-void gmx_tng_prepare_low_prec_writing(gmx_tng_trajectory_t gmx_tng,
- const gmx_mtop_t *mtop,
- const t_inputrec *ir)
+void gmx_tng_prepare_low_prec_writing(gmx_tng_trajectory_t gmx_tng, const gmx_mtop_t* mtop, const t_inputrec* ir)
{
#if GMX_USE_TNG
gmx_tng_add_mtop(gmx_tng, mtop);
add_selection_groups(gmx_tng, mtop);
set_writing_intervals(gmx_tng, TRUE, ir);
- tng_time_per_frame_set(gmx_tng->tng, ir->delta_t * PICO);
+ tng_time_per_frame_set(gmx_tng->tng, ir->delta_t * gmx::c_pico);
gmx_tng->timePerFrameIsSet = true;
gmx_tng_set_compression_precision(gmx_tng, ir->x_compression_precision);
#else
void gmx_fwrite_tng(gmx_tng_trajectory_t gmx_tng,
const gmx_bool bUseLossyCompression,
- gmx_int64_t step,
+ int64_t step,
real elapsedPicoSeconds,
real lambda,
- const rvec *box,
+ const rvec* box,
int nAtoms,
- const rvec *x,
- const rvec *v,
- const rvec *f)
+ const rvec* x,
+ const rvec* v,
+ const rvec* f)
{
#if GMX_USE_TNG
typedef tng_function_status (*write_data_func_pointer)(tng_trajectory_t,
- const gmx_int64_t,
+ const int64_t,
const double,
const real*,
- const gmx_int64_t,
- const gmx_int64_t,
+ const int64_t,
+ const int64_t,
const char*,
const char,
const char);
-#if GMX_DOUBLE
- static write_data_func_pointer write_data = tng_util_generic_with_time_double_write;
-#else
- static write_data_func_pointer write_data = tng_util_generic_with_time_write;
-#endif
- double elapsedSeconds = elapsedPicoSeconds * PICO;
- gmx_int64_t nParticles;
- char compression;
+# if GMX_DOUBLE
+ static write_data_func_pointer write_data = tng_util_generic_with_time_double_write;
+# else
+ static write_data_func_pointer write_data = tng_util_generic_with_time_write;
+# endif
+ double elapsedSeconds = elapsedPicoSeconds * gmx::c_pico;
+ int64_t nParticles;
+ char compression;
if (!gmx_tng)
{
double deltaTime = elapsedSeconds - gmx_tng->lastTime;
std::int64_t deltaStep = step - gmx_tng->lastStep;
- tng_time_per_frame_set(tng, deltaTime / deltaStep );
+ tng_time_per_frame_set(tng, deltaTime / deltaStep);
gmx_tng->timePerFrameIsSet = true;
}
tng_num_particles_get(tng, &nParticles);
- if (nAtoms != (int)nParticles)
+ if (nAtoms != static_cast<int>(nParticles))
{
tng_implicit_num_particles_set(tng, nAtoms);
}
{
GMX_ASSERT(box, "Need a non-NULL box if positions are written");
- if (write_data(tng, step, elapsedSeconds,
- reinterpret_cast<const real *>(x),
- 3, TNG_TRAJ_POSITIONS, "POSITIONS",
+ if (write_data(tng,
+ step,
+ elapsedSeconds,
+ reinterpret_cast<const real*>(x),
+ 3,
+ TNG_TRAJ_POSITIONS,
+ "POSITIONS",
TNG_PARTICLE_BLOCK_DATA,
- compression) != TNG_SUCCESS)
+ compression)
+ != TNG_SUCCESS)
{
gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
}
if (v)
{
- if (write_data(tng, step, elapsedSeconds,
- reinterpret_cast<const real *>(v),
- 3, TNG_TRAJ_VELOCITIES, "VELOCITIES",
+ if (write_data(tng,
+ step,
+ elapsedSeconds,
+ reinterpret_cast<const real*>(v),
+ 3,
+ TNG_TRAJ_VELOCITIES,
+ "VELOCITIES",
TNG_PARTICLE_BLOCK_DATA,
- compression) != TNG_SUCCESS)
+ compression)
+ != TNG_SUCCESS)
{
gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
}
{
/* TNG-MF1 compression only compresses positions and velocities. Use lossless
* compression for forces regardless of output mode */
- if (write_data(tng, step, elapsedSeconds,
- reinterpret_cast<const real *>(f),
- 3, TNG_TRAJ_FORCES, "FORCES",
+ if (write_data(tng,
+ step,
+ elapsedSeconds,
+ reinterpret_cast<const real*>(f),
+ 3,
+ TNG_TRAJ_FORCES,
+ "FORCES",
TNG_PARTICLE_BLOCK_DATA,
- TNG_GZIP_COMPRESSION) != TNG_SUCCESS)
+ TNG_GZIP_COMPRESSION)
+ != TNG_SUCCESS)
{
gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
}
}
- /* TNG-MF1 compression only compresses positions and velocities. Use lossless
- * compression for lambdas and box shape regardless of output mode */
- if (write_data(tng, step, elapsedSeconds,
- reinterpret_cast<const real *>(box),
- 9, TNG_TRAJ_BOX_SHAPE, "BOX SHAPE",
- TNG_NON_PARTICLE_BLOCK_DATA,
- TNG_GZIP_COMPRESSION) != TNG_SUCCESS)
+ if (box)
{
- gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
+ /* TNG-MF1 compression only compresses positions and velocities. Use lossless
+ * compression for box shape regardless of output mode */
+ if (write_data(tng,
+ step,
+ elapsedSeconds,
+ reinterpret_cast<const real*>(box),
+ 9,
+ TNG_TRAJ_BOX_SHAPE,
+ "BOX SHAPE",
+ TNG_NON_PARTICLE_BLOCK_DATA,
+ TNG_GZIP_COMPRESSION)
+ != TNG_SUCCESS)
+ {
+ gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
+ }
}
- if (write_data(tng, step, elapsedSeconds,
- reinterpret_cast<const real *>(&lambda),
- 1, TNG_GMX_LAMBDA, "LAMBDAS",
- TNG_NON_PARTICLE_BLOCK_DATA,
- TNG_GZIP_COMPRESSION) != TNG_SUCCESS)
+ if (lambda >= 0)
{
- gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
+ /* TNG-MF1 compression only compresses positions and velocities. Use lossless
+ * compression for lambda regardless of output mode */
+ if (write_data(tng,
+ step,
+ elapsedSeconds,
+ reinterpret_cast<const real*>(&lambda),
+ 1,
+ TNG_GMX_LAMBDA,
+ "LAMBDAS",
+ TNG_NON_PARTICLE_BLOCK_DATA,
+ TNG_GZIP_COMPRESSION)
+ != TNG_SUCCESS)
+ {
+ gmx_file("Cannot write TNG trajectory frame; maybe you are out of disk space?");
+ }
}
// Update the data in the wrapper
float gmx_tng_get_time_of_final_frame(gmx_tng_trajectory_t gmx_tng)
{
#if GMX_USE_TNG
- gmx_int64_t nFrames;
+ int64_t nFrames;
double time;
float fTime;
tng_trajectory_t tng = gmx_tng->tng;
tng_num_frames_get(tng, &nFrames);
tng_util_time_of_frame_get(tng, nFrames - 1, &time);
- fTime = time / PICO;
+ fTime = time / gmx::c_pico;
return fTime;
#else
GMX_UNUSED_VALUE(gmx_tng);
#endif
}
-void gmx_prepare_tng_writing(const char *filename,
+void gmx_prepare_tng_writing(const char* filename,
char mode,
- gmx_tng_trajectory_t *gmx_tng_input,
- gmx_tng_trajectory_t *gmx_tng_output,
+ gmx_tng_trajectory_t* gmx_tng_input,
+ gmx_tng_trajectory_t* gmx_tng_output,
int nAtoms,
- const gmx_mtop_t *mtop,
- const int *index,
- const char *indexGroupName)
+ const gmx_mtop_t* mtop,
+ gmx::ArrayRef<const int> index,
+ const char* indexGroupName)
{
#if GMX_USE_TNG
- tng_trajectory_t *input = (gmx_tng_input && *gmx_tng_input) ? &(*gmx_tng_input)->tng : nullptr;
+ tng_trajectory_t* input = (gmx_tng_input && *gmx_tng_input) ? &(*gmx_tng_input)->tng : nullptr;
/* FIXME after 5.0: Currently only standard block types are read */
- const int defaultNumIds = 5;
- static gmx_int64_t fallbackIds[defaultNumIds] =
- {
- TNG_TRAJ_BOX_SHAPE, TNG_TRAJ_POSITIONS,
- TNG_TRAJ_VELOCITIES, TNG_TRAJ_FORCES,
- TNG_GMX_LAMBDA
+ const int defaultNumIds = 5;
+ static int64_t fallbackIds[defaultNumIds] = {
+ TNG_TRAJ_BOX_SHAPE, TNG_TRAJ_POSITIONS, TNG_TRAJ_VELOCITIES, TNG_TRAJ_FORCES, TNG_GMX_LAMBDA
};
- static char fallbackNames[defaultNumIds][32] =
- {
- "BOX SHAPE", "POSITIONS", "VELOCITIES",
- "FORCES", "LAMBDAS"
+ static char fallbackNames[defaultNumIds][32] = {
+ "BOX SHAPE", "POSITIONS", "VELOCITIES", "FORCES", "LAMBDAS"
};
- typedef tng_function_status (*set_writing_interval_func_pointer)(tng_trajectory_t,
- const gmx_int64_t,
- const gmx_int64_t,
- const gmx_int64_t,
- const char*,
- const char,
- const char);
-#if GMX_DOUBLE
+ typedef tng_function_status (*set_writing_interval_func_pointer)(
+ tng_trajectory_t, const int64_t, const int64_t, const int64_t, const char*, const char, const char);
+# if GMX_DOUBLE
set_writing_interval_func_pointer set_writing_interval = tng_util_generic_write_interval_double_set;
-#else
+# else
set_writing_interval_func_pointer set_writing_interval = tng_util_generic_write_interval_set;
-#endif
+# endif
gmx_tng_open(filename, mode, gmx_tng_output);
- tng_trajectory_t *output = &(*gmx_tng_output)->tng;
+ tng_trajectory_t* output = &(*gmx_tng_output)->tng;
/* Do we have an input file in TNG format? If so, then there's
more data we can copy over, rather than having to improvise. */
* intervals of positions, box shape and lambdas) of the
* output tng container based on their respective values int
* the input tng container */
- double time, compression_precision;
- gmx_int64_t n_frames_per_frame_set, interval = -1;
+ double time, compression_precision;
+ int64_t n_frames_per_frame_set, interval = -1;
tng_compression_precision_get(*input, &compression_precision);
tng_compression_precision_set(*output, compression_precision);
tng_molecule_system_copy(*input, *output);
tng_time_per_frame_get(*input, &time);
- tng_time_per_frame_set(*output, time);
- // Since we have copied the value from the input TNG we should not change it again
- (*gmx_tng_output)->timePerFrameIsSet = true;
+ /* Only write the time per frame if it was written (and valid). E.g., single
+ * frame files do not usually contain any time per frame information. */
+ if (time >= 0)
+ {
+ (*gmx_tng_input)->timePerFrameIsSet = true;
+ tng_time_per_frame_set(*output, time);
+ // Since we have copied the value from the input TNG we should not change it again
+ (*gmx_tng_output)->timePerFrameIsSet = true;
+ }
tng_num_frames_per_frame_set_get(*input, &n_frames_per_frame_set);
tng_num_frames_per_frame_set_set(*output, n_frames_per_frame_set);
for (int i = 0; i < defaultNumIds; i++)
{
- if (tng_data_get_stride_length(*input, fallbackIds[i], -1, &interval)
- == TNG_SUCCESS)
+ if (tng_data_get_stride_length(*input, fallbackIds[i], -1, &interval) == TNG_SUCCESS)
{
switch (fallbackIds[i])
{
case TNG_TRAJ_POSITIONS:
case TNG_TRAJ_VELOCITIES:
- set_writing_interval(*output, interval, 3, fallbackIds[i],
- fallbackNames[i], TNG_PARTICLE_BLOCK_DATA,
+ set_writing_interval(*output,
+ interval,
+ 3,
+ fallbackIds[i],
+ fallbackNames[i],
+ TNG_PARTICLE_BLOCK_DATA,
compression_type);
break;
case TNG_TRAJ_FORCES:
- set_writing_interval(*output, interval, 3, fallbackIds[i],
- fallbackNames[i], TNG_PARTICLE_BLOCK_DATA,
+ set_writing_interval(*output,
+ interval,
+ 3,
+ fallbackIds[i],
+ fallbackNames[i],
+ TNG_PARTICLE_BLOCK_DATA,
TNG_GZIP_COMPRESSION);
break;
case TNG_TRAJ_BOX_SHAPE:
- set_writing_interval(*output, interval, 9, fallbackIds[i],
- fallbackNames[i], TNG_NON_PARTICLE_BLOCK_DATA,
+ set_writing_interval(*output,
+ interval,
+ 9,
+ fallbackIds[i],
+ fallbackNames[i],
+ TNG_NON_PARTICLE_BLOCK_DATA,
TNG_GZIP_COMPRESSION);
+ (*gmx_tng_output)->boxOutputInterval = interval;
break;
case TNG_GMX_LAMBDA:
- set_writing_interval(*output, interval, 1, fallbackIds[i],
- fallbackNames[i], TNG_NON_PARTICLE_BLOCK_DATA,
+ set_writing_interval(*output,
+ interval,
+ 1,
+ fallbackIds[i],
+ fallbackNames[i],
+ TNG_NON_PARTICLE_BLOCK_DATA,
TNG_GZIP_COMPRESSION);
+ (*gmx_tng_output)->lambdaOutputInterval = interval;
break;
- default:
- continue;
+ default: continue;
}
}
}
-
}
else
{
tng_num_frames_per_frame_set_set(*output, 1);
}
- if (index && nAtoms > 0)
+ if ((!index.empty()) && nAtoms > 0)
{
- gmx_tng_setup_atom_subgroup(*gmx_tng_output, nAtoms, index, indexGroupName);
+ gmx_tng_setup_atom_subgroup(*gmx_tng_output, index, indexGroupName);
}
/* If for some reason there are more requested atoms than there are atoms in the
#endif
}
-void gmx_write_tng_from_trxframe(gmx_tng_trajectory_t gmx_tng_output,
- const t_trxframe *frame,
- int natoms)
+void gmx_write_tng_from_trxframe(gmx_tng_trajectory_t gmx_tng_output, const t_trxframe* frame, int natoms)
{
#if GMX_USE_TNG
if (natoms < 0)
{
natoms = frame->natoms;
}
- gmx_fwrite_tng(gmx_tng_output,
- TRUE,
- frame->step,
- frame->time,
- 0,
- frame->box,
- natoms,
- frame->x,
- frame->v,
- frame->f);
+ gmx_fwrite_tng(
+ gmx_tng_output, TRUE, frame->step, frame->time, 0, frame->box, natoms, frame->x, frame->v, frame->f);
#else
GMX_UNUSED_VALUE(gmx_tng_output);
GMX_UNUSED_VALUE(frame);
{
#if GMX_USE_TNG
-void
-convert_array_to_real_array(void *from,
- real *to,
- const float fact,
- const int nAtoms,
- const int nValues,
- const char datatype)
+void convert_array_to_real_array(void* from,
+ real* to,
+ const float fact,
+ const int nAtoms,
+ const int nValues,
+ const char datatype)
{
- int i, j;
+ int i, j;
const bool useDouble = GMX_DOUBLE;
switch (datatype)
{
for (j = 0; j < nValues; j++)
{
- to[i*nValues+j] = reinterpret_cast<float *>(from)[i*nValues+j] * fact;
+ to[i * nValues + j] = reinterpret_cast<float*>(from)[i * nValues + j] * fact;
}
}
}
{
for (j = 0; j < nValues; j++)
{
- to[i*nValues+j] = reinterpret_cast<float *>(from)[i*nValues+j] * fact;
+ to[i * nValues + j] = reinterpret_cast<float*>(from)[i * nValues + j] * fact;
}
}
}
{
for (j = 0; j < nValues; j++)
{
- to[i*nValues+j] = reinterpret_cast<gmx_int64_t *>(from)[i*nValues+j] * fact;
+ to[i * nValues + j] = reinterpret_cast<int64_t*>(from)[i * nValues + j] * fact;
}
}
break;
{
for (j = 0; j < nValues; j++)
{
- to[i*nValues+j] = reinterpret_cast<double *>(from)[i*nValues+j] * fact;
+ to[i * nValues + j] = reinterpret_cast<double*>(from)[i * nValues + j] * fact;
}
}
}
{
for (j = 0; j < nValues; j++)
{
- to[i*nValues+j] = reinterpret_cast<double *>(from)[i*nValues+j] * fact;
+ to[i * nValues + j] = reinterpret_cast<double*>(from)[i * nValues + j] * fact;
}
}
}
break;
- default:
- gmx_incons("Illegal datatype when converting values to a real array!");
- return;
+ default: gmx_incons("Illegal datatype when converting values to a real array!");
}
- return;
}
real getDistanceScaleFactor(gmx_tng_trajectory_t in)
{
- gmx_int64_t exp = -1;
- real distanceScaleFactor;
+ int64_t exp = -1;
+ real distanceScaleFactor;
// TODO Hopefully, TNG 2.0 will do this kind of thing for us
tng_distance_unit_exponential_get(in->tng, &exp);
// GROMACS expects distances in nm
switch (exp)
{
- case 9:
- distanceScaleFactor = NANO/NANO;
- break;
- case 10:
- distanceScaleFactor = NANO/ANGSTROM;
- break;
- default:
- distanceScaleFactor = pow(10.0, exp + 9.0);
+ case 9: distanceScaleFactor = gmx::c_nano / gmx::c_nano; break;
+ case 10: distanceScaleFactor = gmx::c_nano / gmx::c_angstrom; break;
+ default: distanceScaleFactor = pow(10.0, exp + 9.0);
}
return distanceScaleFactor;
} // namespace
-void gmx_tng_setup_atom_subgroup(gmx_tng_trajectory_t gmx_tng,
- const int nind,
- const int *ind,
- const char *name)
+void gmx_tng_setup_atom_subgroup(gmx_tng_trajectory_t gmx_tng, gmx::ArrayRef<const int> ind, const char* name)
{
#if GMX_USE_TNG
- gmx_int64_t nAtoms, cnt, nMols;
- tng_molecule_t mol, iterMol;
- tng_chain_t chain;
- tng_residue_t res;
- tng_atom_t atom;
- tng_function_status stat;
- tng_trajectory_t tng = gmx_tng->tng;
+ int64_t nAtoms, cnt, nMols;
+ tng_molecule_t mol, iterMol;
+ tng_chain_t chain;
+ tng_residue_t res;
+ tng_atom_t atom;
+ tng_function_status stat;
+ tng_trajectory_t tng = gmx_tng->tng;
tng_num_particles_get(tng, &nAtoms);
- if (nAtoms == nind)
+ if (nAtoms == ind.ssize())
{
return;
}
{
tng_molecule_num_atoms_get(tng, mol, &nAtoms);
tng_molecule_cnt_get(tng, mol, &cnt);
- if (nAtoms == nind)
+ if (nAtoms == ind.ssize())
{
stat = TNG_SUCCESS;
}
tng_molecule_name_set(tng, mol, name);
tng_molecule_chain_add(tng, mol, "", &chain);
- for (int i = 0; i < nind; i++)
+ for (gmx::index i = 0; i < ind.ssize(); i++)
{
- char temp_name[256], temp_type[256];
+ char temp_name[256], temp_type[256];
/* Try to retrieve the residue name of the atom */
stat = tng_residue_name_of_particle_nr_get(tng, ind[i], temp_name, 256);
temp_name[0] = '\0';
}
/* Check if the molecule of the selection already contains this residue */
- if (tng_chain_residue_find(tng, chain, temp_name, -1, &res)
- != TNG_SUCCESS)
+ if (tng_chain_residue_find(tng, chain, temp_name, -1, &res) != TNG_SUCCESS)
{
tng_chain_residue_add(tng, chain, temp_name, &res);
}
* other molecules to 0 */
tng_molecule_cnt_set(tng, mol, 1);
tng_num_molecule_types_get(tng, &nMols);
- for (gmx_int64_t k = 0; k < nMols; k++)
+ for (int64_t k = 0; k < nMols; k++)
{
tng_molecule_of_index_get(tng, k, &iterMol);
if (iterMol == mol)
}
#else
GMX_UNUSED_VALUE(gmx_tng);
- GMX_UNUSED_VALUE(nind);
GMX_UNUSED_VALUE(ind);
GMX_UNUSED_VALUE(name);
#endif
* uncompressing them, then this implemenation should be reconsidered.
* Ideally, gmx trjconv -f a.tng -o b.tng -b 10 -e 20 would be fast
* and lose no information. */
-gmx_bool gmx_read_next_tng_frame(gmx_tng_trajectory_t gmx_tng_input,
- t_trxframe *fr,
- gmx_int64_t *requestedIds,
- int numRequestedIds)
+gmx_bool gmx_read_next_tng_frame(gmx_tng_trajectory_t gmx_tng_input,
+ t_trxframe* fr,
+ int64_t* requestedIds,
+ int numRequestedIds)
{
#if GMX_USE_TNG
- tng_trajectory_t input = gmx_tng_input->tng;
- gmx_bool bOK = TRUE;
- tng_function_status stat;
- gmx_int64_t numberOfAtoms = -1, frameNumber = -1;
- gmx_int64_t nBlocks, blockId, *blockIds = nullptr, codecId;
- char datatype = -1;
- void *values = nullptr;
- double frameTime = -1.0;
- int size, blockDependency;
- double prec;
- const int defaultNumIds = 5;
- static gmx_int64_t fallbackRequestedIds[defaultNumIds] =
- {
- TNG_TRAJ_BOX_SHAPE, TNG_TRAJ_POSITIONS,
- TNG_TRAJ_VELOCITIES, TNG_TRAJ_FORCES,
- TNG_GMX_LAMBDA
+ tng_trajectory_t input = gmx_tng_input->tng;
+ gmx_bool bOK = TRUE;
+ tng_function_status stat;
+ int64_t numberOfAtoms = -1, frameNumber = -1;
+ int64_t nBlocks, blockId, *blockIds = nullptr, codecId;
+ char datatype = -1;
+ void* values = nullptr;
+ double frameTime = -1.0;
+ int size, blockDependency;
+ double prec;
+ const int defaultNumIds = 5;
+ static int64_t fallbackRequestedIds[defaultNumIds] = {
+ TNG_TRAJ_BOX_SHAPE, TNG_TRAJ_POSITIONS, TNG_TRAJ_VELOCITIES, TNG_TRAJ_FORCES, TNG_GMX_LAMBDA
};
- fr->bStep = FALSE;
- fr->bTime = FALSE;
- fr->bLambda = FALSE;
- fr->bAtoms = FALSE;
- fr->bPrec = FALSE;
- fr->bX = FALSE;
- fr->bV = FALSE;
- fr->bF = FALSE;
- fr->bBox = FALSE;
+ fr->bStep = FALSE;
+ fr->bTime = FALSE;
+ fr->bLambda = FALSE;
+ fr->bAtoms = FALSE;
+ fr->bPrec = FALSE;
+ fr->bX = FALSE;
+ fr->bV = FALSE;
+ fr->bF = FALSE;
+ fr->bBox = FALSE;
/* If no specific IDs were requested read all block types that can
* currently be interpreted */
}
fr->natoms = numberOfAtoms;
- bool nextFrameExists = gmx_get_tng_data_block_types_of_next_frame(gmx_tng_input,
- fr->step,
- numRequestedIds,
- requestedIds,
- &frameNumber,
- &nBlocks,
- &blockIds);
- gmx::unique_cptr<gmx_int64_t, gmx::free_wrapper> blockIdsGuard(blockIds);
+ bool nextFrameExists = gmx_get_tng_data_block_types_of_next_frame(
+ gmx_tng_input, fr->step, numRequestedIds, requestedIds, &frameNumber, &nBlocks, &blockIds);
+ gmx::unique_cptr<int64_t, gmx::free_wrapper> blockIdsGuard(blockIds);
if (!nextFrameExists)
{
return FALSE;
return FALSE;
}
- for (gmx_int64_t i = 0; i < nBlocks; i++)
+ for (int64_t i = 0; i < nBlocks; i++)
{
blockId = blockIds[i];
tng_data_block_dependency_get(input, blockId, &blockDependency);
if (blockDependency & TNG_PARTICLE_DEPENDENT)
{
- stat = tng_util_particle_data_next_frame_read(input,
- blockId,
- &values,
- &datatype,
- &frameNumber,
- &frameTime);
+ stat = tng_util_particle_data_next_frame_read(
+ input, blockId, &values, &datatype, &frameNumber, &frameTime);
}
else
{
- stat = tng_util_non_particle_data_next_frame_read(input,
- blockId,
- &values,
- &datatype,
- &frameNumber,
- &frameTime);
+ stat = tng_util_non_particle_data_next_frame_read(
+ input, blockId, &values, &datatype, &frameNumber, &frameTime);
}
if (stat == TNG_CRITICAL)
{
case TNG_TRAJ_BOX_SHAPE:
switch (datatype)
{
- case TNG_INT_DATA:
- size = sizeof(gmx_int64_t);
- break;
- case TNG_FLOAT_DATA:
- size = sizeof(float);
- break;
- case TNG_DOUBLE_DATA:
- size = sizeof(double);
- break;
- default:
- gmx_incons("Illegal datatype of box shape values!");
+ case TNG_INT_DATA: size = sizeof(int64_t); break;
+ case TNG_FLOAT_DATA: size = sizeof(float); break;
+ case TNG_DOUBLE_DATA: size = sizeof(double); break;
+ default: gmx_incons("Illegal datatype of box shape values!");
}
for (int i = 0; i < DIM; i++)
{
- convert_array_to_real_array(reinterpret_cast<char *>(values) + size * i * DIM,
- reinterpret_cast<real *>(fr->box[i]),
+ convert_array_to_real_array(reinterpret_cast<char*>(values) + size * i * DIM,
+ reinterpret_cast<real*>(fr->box[i]),
getDistanceScaleFactor(gmx_tng_input),
1,
DIM,
case TNG_TRAJ_POSITIONS:
srenew(fr->x, fr->natoms);
convert_array_to_real_array(values,
- reinterpret_cast<real *>(fr->x),
+ reinterpret_cast<real*>(fr->x),
getDistanceScaleFactor(gmx_tng_input),
fr->natoms,
DIM,
case TNG_TRAJ_VELOCITIES:
srenew(fr->v, fr->natoms);
convert_array_to_real_array(values,
- (real *) fr->v,
+ reinterpret_cast<real*>(fr->v),
getDistanceScaleFactor(gmx_tng_input),
fr->natoms,
DIM,
case TNG_TRAJ_FORCES:
srenew(fr->f, fr->natoms);
convert_array_to_real_array(values,
- reinterpret_cast<real *>(fr->f),
+ reinterpret_cast<real*>(fr->f),
getDistanceScaleFactor(gmx_tng_input),
fr->natoms,
DIM,
case TNG_GMX_LAMBDA:
switch (datatype)
{
- case TNG_FLOAT_DATA:
- fr->lambda = *(reinterpret_cast<float *>(values));
- break;
- case TNG_DOUBLE_DATA:
- fr->lambda = *(reinterpret_cast<double *>(values));
- break;
- default:
- gmx_incons("Illegal datatype lambda value!");
+ case TNG_FLOAT_DATA: fr->lambda = *(reinterpret_cast<float*>(values)); break;
+ case TNG_DOUBLE_DATA: fr->lambda = *(reinterpret_cast<double*>(values)); break;
+ default: gmx_incons("Illegal datatype lambda value!");
}
fr->bLambda = TRUE;
break;
default:
- gmx_warning("Illegal block type! Currently GROMACS tools can only handle certain data types. Skipping block.");
+ gmx_warning(
+ "Illegal block type! Currently GROMACS tools can only handle certain data "
+ "types. Skipping block.");
}
/* values does not have to be freed before reading next frame. It will
* be reallocated if it is not NULL. */
fr->bStep = TRUE;
// Convert the time to ps
- fr->time = frameTime / PICO;
+ fr->time = frameTime / gmx::c_pico;
fr->bTime = (frameTime > 0);
// TODO This does not leak, but is not exception safe.
#endif
}
-void gmx_print_tng_molecule_system(gmx_tng_trajectory_t gmx_tng_input,
- FILE *stream)
+void gmx_print_tng_molecule_system(gmx_tng_trajectory_t gmx_tng_input, FILE* stream)
{
#if GMX_USE_TNG
- gmx_int64_t nMolecules, nChains, nResidues, nAtoms, nFramesRead;
- gmx_int64_t strideLength, nParticlesRead, nValuesPerFrameRead, *molCntList;
+ int64_t nMolecules, nChains, nResidues, nAtoms, nFramesRead;
+ int64_t strideLength, nParticlesRead, nValuesPerFrameRead, *molCntList;
tng_molecule_t molecule;
tng_chain_t chain;
tng_residue_t residue;
char str[256];
char varNAtoms;
char datatype;
- void *data = nullptr;
+ void* data = nullptr;
std::vector<real> atomCharges;
std::vector<real> atomMasses;
tng_trajectory_t input = gmx_tng_input->tng;
/* Can the number of particles change in the trajectory or is it constant? */
tng_num_particles_variable_get(input, &varNAtoms);
- for (gmx_int64_t i = 0; i < nMolecules; i++)
+ for (int64_t i = 0; i < nMolecules; i++)
{
tng_molecule_of_index_get(input, i, &molecule);
tng_molecule_name_get(input, molecule, str, 256);
if (varNAtoms == TNG_CONSTANT_N_ATOMS)
{
- if ((int)molCntList[i] == 0)
+ if (static_cast<int>(molCntList[i]) == 0)
{
continue;
}
- fprintf(stream, "Molecule: %s, count: %d\n", str, (int)molCntList[i]);
+ fprintf(stream, "Molecule: %s, count: %d\n", str, static_cast<int>(molCntList[i]));
}
else
{
tng_molecule_num_chains_get(input, molecule, &nChains);
if (nChains > 0)
{
- for (gmx_int64_t j = 0; j < nChains; j++)
+ for (int64_t j = 0; j < nChains; j++)
{
tng_molecule_chain_of_index_get(input, molecule, j, &chain);
tng_chain_name_get(input, chain, str, 256);
fprintf(stream, "\tChain: %s\n", str);
tng_chain_num_residues_get(input, chain, &nResidues);
- for (gmx_int64_t k = 0; k < nResidues; k++)
+ for (int64_t k = 0; k < nResidues; k++)
{
tng_chain_residue_of_index_get(input, chain, k, &residue);
tng_residue_name_get(input, residue, str, 256);
fprintf(stream, "\t\tResidue: %s\n", str);
tng_residue_num_atoms_get(input, residue, &nAtoms);
- for (gmx_int64_t l = 0; l < nAtoms; l++)
+ for (int64_t l = 0; l < nAtoms; l++)
{
tng_residue_atom_of_index_get(input, residue, l, &atom);
tng_atom_name_get(input, atom, str, 256);
tng_molecule_num_residues_get(input, molecule, &nResidues);
if (nResidues > 0)
{
- for (gmx_int64_t k = 0; k < nResidues; k++)
+ for (int64_t k = 0; k < nResidues; k++)
{
tng_molecule_residue_of_index_get(input, molecule, k, &residue);
tng_residue_name_get(input, residue, str, 256);
fprintf(stream, "\t\tResidue: %s\n", str);
tng_residue_num_atoms_get(input, residue, &nAtoms);
- for (gmx_int64_t l = 0; l < nAtoms; l++)
+ for (int64_t l = 0; l < nAtoms; l++)
{
tng_residue_atom_of_index_get(input, residue, l, &atom);
tng_atom_name_get(input, atom, str, 256);
else
{
tng_molecule_num_atoms_get(input, molecule, &nAtoms);
- for (gmx_int64_t l = 0; l < nAtoms; l++)
+ for (int64_t l = 0; l < nAtoms; l++)
{
tng_molecule_atom_of_index_get(input, molecule, l, &atom);
tng_atom_name_get(input, atom, str, 256);
}
tng_num_particles_get(input, &nAtoms);
- stat = tng_particle_data_vector_get(input, TNG_TRAJ_PARTIAL_CHARGES, &data, &nFramesRead,
- &strideLength, &nParticlesRead,
- &nValuesPerFrameRead, &datatype);
+ stat = tng_particle_data_vector_get(input,
+ TNG_TRAJ_PARTIAL_CHARGES,
+ &data,
+ &nFramesRead,
+ &strideLength,
+ &nParticlesRead,
+ &nValuesPerFrameRead,
+ &datatype);
if (stat == TNG_SUCCESS)
{
atomCharges.resize(nAtoms);
- convert_array_to_real_array(data,
- atomCharges.data(),
- 1,
- nAtoms,
- 1,
- datatype);
+ convert_array_to_real_array(data, atomCharges.data(), 1, nAtoms, 1, datatype);
fprintf(stream, "Atom Charges (%d):\n", int(nAtoms));
- for (gmx_int64_t i = 0; i < nAtoms; i += 10)
+ for (int64_t i = 0; i < nAtoms; i += 10)
{
fprintf(stream, "Atom Charges [%8d-]=[", int(i));
- for (gmx_int64_t j = 0; (j < 10 && i + j < nAtoms); j++)
+ for (int64_t j = 0; (j < 10 && i + j < nAtoms); j++)
{
fprintf(stream, " %12.5e", atomCharges[i + j]);
}
}
}
- stat = tng_particle_data_vector_get(input, TNG_TRAJ_MASSES, &data, &nFramesRead,
- &strideLength, &nParticlesRead,
- &nValuesPerFrameRead, &datatype);
+ stat = tng_particle_data_vector_get(
+ input, TNG_TRAJ_MASSES, &data, &nFramesRead, &strideLength, &nParticlesRead, &nValuesPerFrameRead, &datatype);
if (stat == TNG_SUCCESS)
{
atomMasses.resize(nAtoms);
- convert_array_to_real_array(data,
- atomMasses.data(),
- 1,
- nAtoms,
- 1,
- datatype);
+ convert_array_to_real_array(data, atomMasses.data(), 1, nAtoms, 1, datatype);
fprintf(stream, "Atom Masses (%d):\n", int(nAtoms));
- for (gmx_int64_t i = 0; i < nAtoms; i += 10)
+ for (int64_t i = 0; i < nAtoms; i += 10)
{
fprintf(stream, "Atom Masses [%8d-]=[", int(i));
- for (gmx_int64_t j = 0; (j < 10 && i + j < nAtoms); j++)
+ for (int64_t j = 0; (j < 10 && i + j < nAtoms); j++)
{
fprintf(stream, " %12.5e", atomMasses[i + j]);
}
gmx_bool gmx_get_tng_data_block_types_of_next_frame(gmx_tng_trajectory_t gmx_tng_input,
int frame,
int nRequestedIds,
- gmx_int64_t *requestedIds,
- gmx_int64_t *nextFrame,
- gmx_int64_t *nBlocks,
- gmx_int64_t **blockIds)
+ int64_t* requestedIds,
+ int64_t* nextFrame,
+ int64_t* nBlocks,
+ int64_t** blockIds)
{
#if GMX_USE_TNG
tng_function_status stat;
tng_trajectory_t input = gmx_tng_input->tng;
- stat = tng_util_trajectory_next_frame_present_data_blocks_find(input, frame,
- nRequestedIds, requestedIds,
- nextFrame,
- nBlocks, blockIds);
+ stat = tng_util_trajectory_next_frame_present_data_blocks_find(
+ input, frame, nRequestedIds, requestedIds, nextFrame, nBlocks, blockIds);
if (stat == TNG_CRITICAL)
{
}
gmx_bool gmx_get_tng_data_next_frame_of_block_type(gmx_tng_trajectory_t gmx_tng_input,
- gmx_int64_t blockId,
- real **values,
- gmx_int64_t *frameNumber,
- double *frameTime,
- gmx_int64_t *nValuesPerFrame,
- gmx_int64_t *nAtoms,
- real *prec,
- char *name,
+ int64_t blockId,
+ real** values,
+ int64_t* frameNumber,
+ double* frameTime,
+ int64_t* nValuesPerFrame,
+ int64_t* nAtoms,
+ real* prec,
+ char* name,
int maxLen,
- gmx_bool *bOK)
+ gmx_bool* bOK)
{
#if GMX_USE_TNG
tng_function_status stat;
char datatype = -1;
- gmx_int64_t codecId;
+ int64_t codecId;
int blockDependency;
- void *data = nullptr;
+ void* data = nullptr;
double localPrec;
tng_trajectory_t input = gmx_tng_input->tng;
if (blockDependency & TNG_PARTICLE_DEPENDENT)
{
tng_num_particles_get(input, nAtoms);
- stat = tng_util_particle_data_next_frame_read(input,
- blockId,
- &data,
- &datatype,
- frameNumber,
- frameTime);
+ stat = tng_util_particle_data_next_frame_read(
+ input, blockId, &data, &datatype, frameNumber, frameTime);
}
else
{
*nAtoms = 1; /* There are not actually any atoms, but it is used for
allocating memory */
- stat = tng_util_non_particle_data_next_frame_read(input,
- blockId,
- &data,
- &datatype,
- frameNumber,
- frameTime);
+ stat = tng_util_non_particle_data_next_frame_read(
+ input, blockId, &data, &datatype, frameNumber, frameTime);
}
if (stat == TNG_CRITICAL)
{
gmx_file("Cannot read next frame of TNG file");
}
srenew(*values, sizeof(real) * *nValuesPerFrame * *nAtoms);
- convert_array_to_real_array(data,
- *values,
- getDistanceScaleFactor(gmx_tng_input),
- *nAtoms,
- *nValuesPerFrame,
- datatype);
+ convert_array_to_real_array(
+ data, *values, getDistanceScaleFactor(gmx_tng_input), *nAtoms, *nValuesPerFrame, datatype);
tng_util_frame_current_compression_get(input, blockId, &codecId, &localPrec);
return FALSE;
#endif
}
+
+int gmx_tng_get_box_output_interval(gmx_tng_trajectory_t gmx_tng)
+{
+#if GMX_USE_TNG
+ return gmx_tng->boxOutputInterval;
+#else
+ GMX_UNUSED_VALUE(gmx_tng);
+ return -1;
+#endif
+}
+
+int gmx_tng_get_lambda_output_interval(gmx_tng_trajectory_t gmx_tng)
+{
+#if GMX_USE_TNG
+ return gmx_tng->lambdaOutputInterval;
+#else
+ GMX_UNUSED_VALUE(gmx_tng);
+ return -1;
+#endif
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff