for (i = 0; i < natoms; i++)
{
ii = dd->gatindex[i];
- mtopGetAtomAndResidueName(mtop, ii, &molb, &atomname, &resnr, &resname, NULL);
+ mtopGetAtomAndResidueName(mtop, ii, &molb, &atomname, &resnr, &resname, nullptr);
if (i < dd->comm->nat[ddnatZONE])
{
c = 0;
snew(pairleg[i], 30);
j = fa[3*i+1];
k = fa[3*i+2];
- mtopGetAtomAndResidueName(&mtop, j, &molb, &anm_j, &resnr_j, &resnm_j, NULL);
- mtopGetAtomAndResidueName(&mtop, k, &molb, &anm_k, &resnr_k, &resnm_k, NULL);
+ mtopGetAtomAndResidueName(&mtop, j, &molb, &anm_j, &resnr_j, &resnm_j, nullptr);
+ mtopGetAtomAndResidueName(&mtop, k, &molb, &anm_k, &resnr_k, &resnm_k, nullptr);
sprintf(pairleg[i], "%d %s %d %s (%d)",
resnr_j, anm_j, resnr_k, anm_k,
ip[fa[3*i]].disres.label);
{
ii = i;
}
- mtopGetAtomAndResidueName(mtop, ii, &molb, &anm, &resnr, &resnm, NULL);
+ mtopGetAtomAndResidueName(mtop, ii, &molb, &anm, &resnr, &resnm, nullptr);
gmx_fprintf_pdb_atomline(out, epdbATOM, ii+1, anm, ' ', resnm, ' ', resnr, ' ',
10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], 1.0, 0.0, "");
}
#include "gromacs/math/vec.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/topology/block.h"
+#include "gromacs/topology/mtop_lookup.h"
#include "gromacs/topology/topology.h"
#include "gromacs/utility/gmxassert.h"
void
-gmx_calc_cog(const t_topology * /* top */, rvec x[], int nrefat, const int index[], rvec xout)
+gmx_calc_cog(const gmx_mtop_t * /* top */, rvec x[], int nrefat, const int index[], rvec xout)
{
int m, ai;
* mass are calculated, and hence a topology with masses is required.
*/
void
-gmx_calc_com(const t_topology *top, rvec x[], int nrefat, const int index[], rvec xout)
+gmx_calc_com(const gmx_mtop_t *top, rvec x[], int nrefat, const int index[], rvec xout)
{
- int m, j, ai;
- real mass, mtot;
-
- GMX_RELEASE_ASSERT(top != nullptr && top->atoms.haveMass,
+ GMX_RELEASE_ASSERT(gmx_mtop_has_masses(top),
"No masses available while mass weighting was requested");
clear_rvec(xout);
- mtot = 0;
- for (m = 0; m < nrefat; ++m)
+ real mtot = 0;
+ int molb = 0;
+ for (int m = 0; m < nrefat; ++m)
{
- ai = index[m];
- mass = top->atoms.atom[ai].m;
- for (j = 0; j < DIM; ++j)
+ const int ai = index[m];
+ const real mass = mtopGetAtomMass(top, ai, &molb);
+ for (int j = 0; j < DIM; ++j)
{
xout[j] += mass * x[ai][j];
}
* \param[out] fout Force on the COG position for the indexed atoms.
*/
void
-gmx_calc_cog_f(const t_topology *top, rvec f[], int nrefat, const int index[], rvec fout)
+gmx_calc_cog_f(const gmx_mtop_t *top, rvec f[], int nrefat, const int index[], rvec fout)
{
- int m, j, ai;
- real mass, mtot;
-
- GMX_RELEASE_ASSERT(top != nullptr && top->atoms.haveMass,
+ GMX_RELEASE_ASSERT(gmx_mtop_has_masses(top),
"No masses available while mass weighting was requested");
clear_rvec(fout);
- mtot = 0;
- for (m = 0; m < nrefat; ++m)
+ real mtot = 0;
+ int molb = 0;
+ for (int m = 0; m < nrefat; ++m)
{
- ai = index[m];
- mass = top->atoms.atom[ai].m;
- for (j = 0; j < DIM; ++j)
+ const int ai = index[m];
+ const real mass = mtopGetAtomMass(top, ai, &molb);
+ for (int j = 0; j < DIM; ++j)
{
fout[j] += f[ai][j] / mass;
}
}
void
-gmx_calc_com_f(const t_topology * /* top */, rvec f[], int nrefat, const int index[], rvec fout)
+gmx_calc_com_f(const gmx_mtop_t * /* top */, rvec f[], int nrefat, const int index[], rvec fout)
{
clear_rvec(fout);
for (int m = 0; m < nrefat; ++m)
* Other parameters are passed unmodified to these functions.
*/
void
-gmx_calc_comg(const t_topology *top, rvec x[], int nrefat, const int index[],
+gmx_calc_comg(const gmx_mtop_t *top, rvec x[], int nrefat, const int index[],
bool bMass, rvec xout)
{
if (bMass)
* Other parameters are passed unmodified to these functions.
*/
void
-gmx_calc_comg_f(const t_topology *top, rvec f[], int nrefat, const int index[],
+gmx_calc_comg_f(const gmx_mtop_t *top, rvec f[], int nrefat, const int index[],
bool bMass, rvec fout)
{
if (bMass)
* Works exactly as gmx_calc_com_pbc(), but calculates the center of geometry.
*/
void
-gmx_calc_cog_pbc(const t_topology *top, rvec x[], const t_pbc *pbc,
+gmx_calc_cog_pbc(const gmx_mtop_t *top, rvec x[], const t_pbc *pbc,
int nrefat, const int index[], rvec xout)
{
const real tol = 1e-4;
* Modified from src/tools/gmx_sorient.c in Gromacs distribution.
*/
void
-gmx_calc_com_pbc(const t_topology *top, rvec x[], const t_pbc *pbc,
+gmx_calc_com_pbc(const gmx_mtop_t *top, rvec x[], const t_pbc *pbc,
int nrefat, const int index[], rvec xout)
{
- const real tol = 1e-4;
- bool bChanged;
- int m, j, ai, iter;
- real mass, mtot;
- rvec dx, xtest;
-
- GMX_RELEASE_ASSERT(top != nullptr && top->atoms.haveMass,
+ GMX_RELEASE_ASSERT(gmx_mtop_has_masses(top),
"No masses available while mass weighting was requested");
/* First simple calculation */
clear_rvec(xout);
- mtot = 0;
- for (m = 0; m < nrefat; ++m)
+ real mtot = 0;
+ int molb = 0;
+ for (int m = 0; m < nrefat; ++m)
{
- ai = index[m];
- mass = top->atoms.atom[ai].m;
- for (j = 0; j < DIM; ++j)
+ const int ai = index[m];
+ const real mass = mtopGetAtomMass(top, ai, &molb);
+ for (int j = 0; j < DIM; ++j)
{
xout[j] += mass * x[ai][j];
}
/* Now check if any atom is more than half the box from the COM */
if (pbc)
{
- iter = 0;
+ const real tol = 1e-4;
+ bool bChanged;
do
{
bChanged = false;
- for (m = 0; m < nrefat; ++m)
+ molb = 0;
+ for (int m = 0; m < nrefat; ++m)
{
- ai = index[m];
- mass = top->atoms.atom[ai].m / mtot;
+ rvec dx, xtest;
+ const int ai = index[m];
+ const real mass = mtopGetAtomMass(top, ai, &molb) / mtot;
pbc_dx(pbc, x[ai], xout, dx);
rvec_add(xout, dx, xtest);
- for (j = 0; j < DIM; ++j)
+ for (int j = 0; j < DIM; ++j)
{
if (fabs(xtest[j] - x[ai][j]) > tol)
{
}
}
}
- iter++;
}
while (bChanged);
}
* Other parameters are passed unmodified to these functions.
*/
void
-gmx_calc_comg_pbc(const t_topology *top, rvec x[], const t_pbc *pbc,
+gmx_calc_comg_pbc(const gmx_mtop_t *top, rvec x[], const t_pbc *pbc,
int nrefat, const int index[], bool bMass, rvec xout)
{
if (bMass)
void
-gmx_calc_cog_block(const t_topology * /* top */, rvec x[], const t_block *block, const int index[],
+gmx_calc_cog_block(const gmx_mtop_t * /* top */, rvec x[], const t_block *block, const int index[],
rvec xout[])
{
int b, i, ai;
* mass are calculated, and hence a topology with masses is required.
*/
void
-gmx_calc_com_block(const t_topology *top, rvec x[], const t_block *block, const int index[],
+gmx_calc_com_block(const gmx_mtop_t *top, rvec x[], const t_block *block, const int index[],
rvec xout[])
{
- int b, i, ai, d;
- rvec xb;
- real mass, mtot;
-
- GMX_RELEASE_ASSERT(top != nullptr && top->atoms.haveMass,
+ GMX_RELEASE_ASSERT(gmx_mtop_has_masses(top),
"No masses available while mass weighting was requested");
- for (b = 0; b < block->nr; ++b)
+ int molb = 0;
+ for (int b = 0; b < block->nr; ++b)
{
+ rvec xb;
clear_rvec(xb);
- mtot = 0;
- for (i = block->index[b]; i < block->index[b+1]; ++i)
+ real mtot = 0;
+ for (int i = block->index[b]; i < block->index[b+1]; ++i)
{
- ai = index[i];
- mass = top->atoms.atom[ai].m;
- for (d = 0; d < DIM; ++d)
+ const int ai = index[i];
+ const real mass = mtopGetAtomMass(top, ai, &molb);
+ for (int d = 0; d < DIM; ++d)
{
xb[d] += mass * x[ai][d];
}
* \param[out] fout \p block->nr Forces on COG positions.
*/
void
-gmx_calc_cog_f_block(const t_topology *top, rvec f[], const t_block *block, const int index[],
+gmx_calc_cog_f_block(const gmx_mtop_t *top, rvec f[], const t_block *block, const int index[],
rvec fout[])
{
- int b, i, ai, d;
- rvec fb;
- real mass, mtot;
-
- GMX_RELEASE_ASSERT(top != nullptr && top->atoms.haveMass,
+ GMX_RELEASE_ASSERT(gmx_mtop_has_masses(top),
"No masses available while mass weighting was requested");
- for (b = 0; b < block->nr; ++b)
+ int molb = 0;
+ for (int b = 0; b < block->nr; ++b)
{
+ rvec fb;
clear_rvec(fb);
- mtot = 0;
- for (i = block->index[b]; i < block->index[b+1]; ++i)
+ real mtot = 0;
+ for (int i = block->index[b]; i < block->index[b+1]; ++i)
{
- ai = index[i];
- mass = top->atoms.atom[ai].m;
- for (d = 0; d < DIM; ++d)
+ const int ai = index[i];
+ const real mass = mtopGetAtomMass(top, ai, &molb);
+ for (int d = 0; d < DIM; ++d)
{
fb[d] += f[ai][d] / mass;
}
}
void
-gmx_calc_com_f_block(const t_topology * /* top */, rvec f[], const t_block *block, const int index[],
+gmx_calc_com_f_block(const gmx_mtop_t * /* top */, rvec f[], const t_block *block, const int index[],
rvec fout[])
{
for (int b = 0; b < block->nr; ++b)
* Other parameters are passed unmodified to these functions.
*/
void
-gmx_calc_comg_block(const t_topology *top, rvec x[], const t_block *block, const int index[],
+gmx_calc_comg_block(const gmx_mtop_t *top, rvec x[], const t_block *block, const int index[],
bool bMass, rvec xout[])
{
if (bMass)
* Other parameters are passed unmodified to these functions.
*/
void
-gmx_calc_comg_f_block(const t_topology *top, rvec f[], const t_block *block, const int index[],
+gmx_calc_comg_f_block(const gmx_mtop_t *top, rvec f[], const t_block *block, const int index[],
bool bMass, rvec fout[])
{
if (bMass)
* crashes.
*/
void
-gmx_calc_comg_blocka(const t_topology *top, rvec x[], const t_blocka *block,
+gmx_calc_comg_blocka(const gmx_mtop_t *top, rvec x[], const t_blocka *block,
bool bMass, rvec xout[])
{
/* TODO: It would probably be better to do this without the type cast */
* crashes.
*/
void
-gmx_calc_comg_f_blocka(const t_topology *top, rvec f[], const t_blocka *block,
+gmx_calc_comg_f_blocka(const gmx_mtop_t *top, rvec f[], const t_blocka *block,
bool bMass, rvec fout[])
{
/* TODO: It would probably be better to do this without the type cast */
#include "gromacs/math/vectypes.h"
+struct gmx_mtop_t;
struct t_block;
struct t_blocka;
struct t_pbc;
-struct t_topology;
/*! \brief
* Calculate a single center of geometry.
* \param[out] xout COG position for the indexed atoms.
*/
void
-gmx_calc_cog(const t_topology *top, rvec x[], int nrefat, const int index[], rvec xout);
+gmx_calc_cog(const gmx_mtop_t *top, rvec x[], int nrefat, const int index[], rvec xout);
/** Calculate a single center of mass. */
void
-gmx_calc_com(const t_topology *top, rvec x[], int nrefat, const int index[], rvec xout);
+gmx_calc_com(const gmx_mtop_t *top, rvec x[], int nrefat, const int index[], rvec xout);
/** Calculate force on a single center of geometry. */
void
-gmx_calc_cog_f(const t_topology *top, rvec f[], int nrefat, const int index[], rvec fout);
+gmx_calc_cog_f(const gmx_mtop_t *top, rvec f[], int nrefat, const int index[], rvec fout);
/*! \brief
* Calculate force on a single center of mass.
*
* \param[out] fout Force on the COM position for the indexed atoms.
*/
void
-gmx_calc_com_f(const t_topology *top, rvec f[], int nrefat, const int index[], rvec fout);
+gmx_calc_com_f(const gmx_mtop_t *top, rvec f[], int nrefat, const int index[], rvec fout);
/** Calculate a single center of mass/geometry. */
void
-gmx_calc_comg(const t_topology *top, rvec x[], int nrefat, const int index[],
+gmx_calc_comg(const gmx_mtop_t *top, rvec x[], int nrefat, const int index[],
bool bMass, rvec xout);
/** Calculate force on a single center of mass/geometry. */
void
-gmx_calc_comg_f(const t_topology *top, rvec f[], int nrefat, const int index[],
+gmx_calc_comg_f(const gmx_mtop_t *top, rvec f[], int nrefat, const int index[],
bool bMass, rvec fout);
/** Calculate a single center of geometry iteratively, taking PBC into account. */
void
-gmx_calc_cog_pbc(const t_topology *top, rvec x[], const t_pbc *pbc,
+gmx_calc_cog_pbc(const gmx_mtop_t *top, rvec x[], const t_pbc *pbc,
int nrefat, const int index[], rvec xout);
/** Calculate a single center of mass iteratively, taking PBC into account. */
void
-gmx_calc_com_pbc(const t_topology *top, rvec x[], const t_pbc *pbc,
+gmx_calc_com_pbc(const gmx_mtop_t *top, rvec x[], const t_pbc *pbc,
int nrefat, const int index[], rvec xout);
/** Calculate a single center of mass/geometry iteratively with PBC. */
void
-gmx_calc_comg_pbc(const t_topology *top, rvec x[], const t_pbc *pbc,
+gmx_calc_comg_pbc(const gmx_mtop_t *top, rvec x[], const t_pbc *pbc,
int nrefat, const int index[], bool bMass, rvec xout);
/*! \brief
* \param[out] xout \p block->nr COG positions.
*/
void
-gmx_calc_cog_block(const t_topology *top, rvec x[], const t_block *block,
+gmx_calc_cog_block(const gmx_mtop_t *top, rvec x[], const t_block *block,
const int index[], rvec xout[]);
/** Calculate centers of mass for a blocked index. */
void
-gmx_calc_com_block(const t_topology *top, rvec x[], const t_block *block,
+gmx_calc_com_block(const gmx_mtop_t *top, rvec x[], const t_block *block,
const int index[], rvec xout[]);
/** Calculate forces on centers of geometry for a blocked index. */
void
-gmx_calc_cog_f_block(const t_topology *top, rvec f[], const t_block *block,
+gmx_calc_cog_f_block(const gmx_mtop_t *top, rvec f[], const t_block *block,
const int index[], rvec fout[]);
/*! \brief
* Calculate forces on centers of mass for a blocked index.
* \param[out] fout \p block->nr Forces on COM positions.
*/
void
-gmx_calc_com_f_block(const t_topology *top, rvec f[], const t_block *block,
+gmx_calc_com_f_block(const gmx_mtop_t *top, rvec f[], const t_block *block,
const int index[], rvec fout[]);
/** Calculate centers of mass/geometry for a blocked index. */
void
-gmx_calc_comg_block(const t_topology *top, rvec x[], const t_block *block,
+gmx_calc_comg_block(const gmx_mtop_t *top, rvec x[], const t_block *block,
const int index[], bool bMass, rvec xout[]);
/** Calculate forces on centers of mass/geometry for a blocked index. */
void
-gmx_calc_comg_f_block(const t_topology *top, rvec f[], const t_block *block,
+gmx_calc_comg_f_block(const gmx_mtop_t *top, rvec f[], const t_block *block,
const int index[], bool bMass, rvec fout[]);
/** Calculate centers of mass/geometry for a set of blocks; */
void
-gmx_calc_comg_blocka(const t_topology *top, rvec x[], const t_blocka *block,
+gmx_calc_comg_blocka(const gmx_mtop_t *top, rvec x[], const t_blocka *block,
bool bMass, rvec xout[]);
/** Calculate forces on centers of mass/geometry for a set of blocks; */
void
-gmx_calc_comg_f_blocka(const t_topology *top, rvec x[], const t_blocka *block,
+gmx_calc_comg_f_blocka(const gmx_mtop_t *top, rvec x[], const t_blocka *block,
bool bMass, rvec xout[]);
#endif
#include "gromacs/topology/block.h"
#include "gromacs/topology/index.h"
+#include "gromacs/topology/mtop_lookup.h"
#include "gromacs/topology/mtop_util.h"
#include "gromacs/topology/topology.h"
#include "gromacs/utility/exceptions.h"
* \ingroup module_selection
*/
static bool
-next_group_index(int atomIndex, t_topology *top, e_index_t type, int *id)
+next_group_index(int atomIndex, const gmx_mtop_t *top, e_index_t type, int *id)
{
int prev = *id;
switch (type)
*id = atomIndex;
break;
case INDEX_RES:
- *id = top->atoms.atom[atomIndex].resind;
+ {
+ int resind, molb = 0;
+ mtopGetAtomAndResidueName(top, atomIndex, &molb, nullptr, nullptr, nullptr, &resind);
+ *id = resind;
break;
+ }
case INDEX_MOL:
if (*id >= 0 && top->mols.index[*id] > atomIndex)
{
* \p g should be sorted.
*/
void
-gmx_ana_index_make_block(t_blocka *t, t_topology *top, gmx_ana_index_t *g,
+gmx_ana_index_make_block(t_blocka *t, const gmx_mtop_t *top, gmx_ana_index_t *g,
e_index_t type, bool bComplete)
{
if (type == INDEX_UNKNOWN)
t->nra = 0;
/* We may allocate some extra memory here because we don't know in
* advance how much will be needed. */
- if (t->nalloc_a < top->atoms.nr)
+ if (t->nalloc_a < top->natoms)
{
- srenew(t->a, top->atoms.nr);
- t->nalloc_a = top->atoms.nr;
+ srenew(t->a, top->natoms);
+ t->nalloc_a = top->natoms;
}
}
else
}
/* Clear counters */
t->nr = 0;
- int j = 0; /* j is used by residue completion for the first atom not stored */
- int id = -1;
+ int id = -1;
+ int molb = 0;
for (int i = 0; i < g->isize; ++i)
{
+ const int ai = g->index[i];
/* Find the ID number of the atom/residue/molecule corresponding to
* the atom. */
- if (next_group_index(g->index[i], top, type, &id))
+ if (next_group_index(ai, top, type, &id))
{
/* If this is the first atom in a new block, initialize the block. */
if (bComplete)
switch (type)
{
case INDEX_RES:
- while (top->atoms.atom[j].resind != id)
+ {
+ int molnr, atnr_mol;
+ mtopGetMolblockIndex(top, ai, &molb, &molnr, &atnr_mol);
+ const t_atoms &mol_atoms = top->moltype[top->molblock[molb].type].atoms;
+ int last_atom = atnr_mol + 1;
+ while (last_atom < mol_atoms.nr
+ && mol_atoms.atom[last_atom].resind == id)
+ {
+ ++last_atom;
+ }
+ int first_atom = atnr_mol - 1;
+ while (first_atom >= 0
+ && mol_atoms.atom[first_atom].resind == id)
{
- ++j;
+ --first_atom;
}
- while (j < top->atoms.nr && top->atoms.atom[j].resind == id)
+ int first_mol_atom = top->molblock[molb].globalAtomStart;
+ first_mol_atom += molnr*top->molblock[molb].natoms_mol;
+ first_atom = first_mol_atom + first_atom + 1;
+ last_atom = first_mol_atom + last_atom - 1;
+ for (int j = first_atom; j <= last_atom; ++j)
{
t->a[t->nra++] = j;
- ++j;
}
break;
-
+ }
case INDEX_MOL:
- for (j = top->mols.index[id]; j < top->mols.index[id+1]; ++j)
+ for (int j = top->mols.index[id]; j < top->mols.index[id+1]; ++j)
{
t->a[t->nra++] = j;
}
* The atoms in \p g are assumed to be sorted.
*/
bool
-gmx_ana_index_has_full_blocks(gmx_ana_index_t *g, t_block *b)
+gmx_ana_index_has_full_blocks(const gmx_ana_index_t *g, const t_block *b)
{
int i, j, bi;
return true;
}
+/*!
+ * \brief Returns if an atom is at a residue boundary.
+ *
+ * \param[in] top Topology data.
+ * \param[in] a Atom index to check, should be -1 <= \p a < top->natoms.
+ * \param[in,out] molb The molecule block of atom a
+ * \returns true if atoms \p a and \p a + 1 are in different residues, false otherwise.
+ */
+static bool is_at_residue_boundary(const gmx_mtop_t *top, int a, int *molb)
+{
+ if (a == -1 || a + 1 == top->natoms)
+ {
+ return true;
+ }
+ int resindA;
+ mtopGetAtomAndResidueName(top, a, molb,
+ nullptr, nullptr, nullptr, &resindA);
+ int resindAPlusOne;
+ mtopGetAtomAndResidueName(top, a + 1, molb,
+ nullptr, nullptr, nullptr, &resindAPlusOne);
+ return resindAPlusOne != resindA;
+}
+
/*!
* \param[in] g Index group to check.
* \param[in] type Block data to check against.
*/
bool
gmx_ana_index_has_complete_elems(gmx_ana_index_t *g, e_index_t type,
- t_topology *top)
+ const gmx_mtop_t *top)
{
+ if (g->isize == 0)
+ {
+ return true;
+ }
+
// TODO: Consider whether unsorted groups need to be supported better.
switch (type)
{
case INDEX_RES:
{
- int i, ai;
- int id, prev;
-
- prev = -1;
- for (i = 0; i < g->isize; ++i)
+ int molb = 0;
+ int aPrev = -1;
+ for (int i = 0; i < g->isize; ++i)
{
- ai = g->index[i];
- id = top->atoms.atom[ai].resind;
- if (id != prev)
+ const int a = g->index[i];
+ // Check if a is consecutive or on a residue boundary
+ if (a != aPrev + 1)
{
- if (ai > 0 && top->atoms.atom[ai-1].resind == id)
+ if (!is_at_residue_boundary(top, aPrev, &molb))
{
return false;
}
- if (i > 0 && g->index[i-1] < top->atoms.nr - 1
- && top->atoms.atom[g->index[i-1]+1].resind == prev)
+ if (!is_at_residue_boundary(top, a - 1, &molb))
{
return false;
}
}
- prev = id;
+ aPrev = a;
}
- if (g->index[i-1] < top->atoms.nr - 1
- && top->atoms.atom[g->index[i-1]+1].resind == prev)
+ GMX_ASSERT(g->isize > 0, "We return above when isize=0");
+ const int a = g->index[g->isize - 1];
+ if (!is_at_residue_boundary(top, a, &molb))
{
return false;
}
*/
void
gmx_ana_indexmap_init(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
- t_topology *top, e_index_t type)
+ const gmx_mtop_t *top, e_index_t type)
{
m->type = type;
gmx_ana_index_make_block(&m->b, top, g, type, false);
}
int
-gmx_ana_indexmap_init_orgid_group(gmx_ana_indexmap_t *m, t_topology *top,
+gmx_ana_indexmap_init_orgid_group(gmx_ana_indexmap_t *m, const gmx_mtop_t *top,
e_index_t type)
{
GMX_RELEASE_ASSERT(m->bStatic,
}
struct gmx_mtop_t;
-struct t_topology;
/** Stores a set of index groups. */
struct gmx_ana_indexgrps_t;
/*@{*/
/** Partition a group based on topology information. */
void
-gmx_ana_index_make_block(t_blocka *t, t_topology *top, gmx_ana_index_t *g,
+gmx_ana_index_make_block(t_blocka *t, const gmx_mtop_t *top, gmx_ana_index_t *g,
e_index_t type, bool bComplete);
/** Checks whether a group consists of full blocks. */
bool
-gmx_ana_index_has_full_blocks(gmx_ana_index_t *g, t_block *b);
+gmx_ana_index_has_full_blocks(const gmx_ana_index_t *g, const t_block *b);
/** Checks whether a group consists of full blocks. */
bool
gmx_ana_index_has_full_ablocks(gmx_ana_index_t *g, t_blocka *b);
/** Checks whether a group consists of full residues/molecules. */
bool
-gmx_ana_index_has_complete_elems(gmx_ana_index_t *g, e_index_t type, t_topology *top);
+gmx_ana_index_has_complete_elems(gmx_ana_index_t *g, e_index_t type,
+ const gmx_mtop_t *top);
/** Initializes an empty index group mapping. */
void
/** Initializes an index group mapping. */
void
gmx_ana_indexmap_init(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
- t_topology *top, e_index_t type);
+ const gmx_mtop_t *top, e_index_t type);
/*! \brief
* Initializes `orgid` entries based on topology grouping.
*
* Strong exception safety guarantee.
*/
int
-gmx_ana_indexmap_init_orgid_group(gmx_ana_indexmap_t *m, t_topology *top,
+gmx_ana_indexmap_init_orgid_group(gmx_ana_indexmap_t *m, const gmx_mtop_t *top,
e_index_t type);
/** Sets an index group mapping to be static. */
void
* Can be NULL if none of the calculations require topology data or if
* setTopology() has not been called.
*/
- t_topology *top_;
+ const gmx_mtop_t *top_;
//! Pointer to the first data structure.
gmx_ana_poscalc_t *first_;
//! Pointer to the last data structure.
}
void
-PositionCalculationCollection::setTopology(t_topology *top)
+PositionCalculationCollection::setTopology(const gmx_mtop_t *top)
{
impl_->top_ = top;
}
static void
set_poscalc_maxindex(gmx_ana_poscalc_t *pc, gmx_ana_index_t *g, bool bBase)
{
- t_topology *top = pc->coll->top_;
+ const gmx_mtop_t *top = pc->coll->top_;
gmx_ana_index_make_block(&pc->b, top, g, pc->itype, pc->flags & POS_COMPLWHOLE);
/* Set the type to POS_ATOM if the calculation in fact is such. */
if (pc->b.nr == pc->b.nra)
}
}
gmx::ConstArrayRef<int> index = pc->coll->getFrameIndices(pc->b.nra, pc->b.a);
- const t_topology *top = pc->coll->top_;
+ const gmx_mtop_t *top = pc->coll->top_;
const bool bMass = pc->flags & POS_MASS;
switch (pc->type)
{
struct gmx_ana_index_t;
struct gmx_ana_pos_t;
+struct gmx_mtop_t;
struct t_pbc;
-struct t_topology;
struct t_trxframe;
namespace gmx
*
* Does not throw.
*/
- void setTopology(t_topology *top);
+ void setTopology(const gmx_mtop_t *top);
/*! \brief
* Prints information about calculations.
*
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2009,2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016, 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.
int
-Selection::initOriginalIdsToGroup(t_topology *top, e_index_t type)
+Selection::initOriginalIdsToGroup(const gmx_mtop_t *top, e_index_t type)
{
try
{
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2009,2010,2011,2012,2013,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016, 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 "gromacs/utility/classhelpers.h"
#include "gromacs/utility/gmxassert.h"
+struct gmx_mtop_t;
struct t_topology;
namespace gmx
* \see setOriginalId()
* \see SelectionPosition::mappedId()
*/
- int initOriginalIdsToGroup(t_topology *top, e_index_t type);
+ int initOriginalIdsToGroup(const gmx_mtop_t *top, e_index_t type);
/*! \brief
* Prints out one-line description of the selection.
{
snew(sc->top, 1);
*sc->top = gmx_mtop_t_to_t_topology(top, false);
- sc->pcc.setTopology(sc->top);
checkTopologyProperties(sc->top, requiredTopologyProperties());
}
else
{
checkTopologyProperties(nullptr, requiredTopologyProperties());
}
+ sc->pcc.setTopology(top);
}
topManager_.initAtoms(15);
topManager_.initUniformResidues(3);
- t_topology *top = topManager_.topology();
+ gmx_mtop_t *top = topManager_.topology();
setGroup(group1);
EXPECT_TRUE(gmx_ana_index_has_complete_elems(&g_, INDEX_RES, top));
const int group1[] = { 3, 4, 5, 6, 7, 8, 12, 13 };
const int group2[] = { 3, 4, 5, 6, 7, 12, 13, 14 };
const int group3[] = { 4, 5, 6, 7, 8, 12, 13, 14 };
+ const int group4[] = { 3, 4, 5, 6, 8, 12, 13, 14 };
topManager_.initAtoms(18);
topManager_.initUniformResidues(3);
- t_topology *top = topManager_.topology();
+ gmx_mtop_t *top = topManager_.topology();
setGroup(group1);
EXPECT_FALSE(gmx_ana_index_has_complete_elems(&g_, INDEX_RES, top));
setGroup(group3);
EXPECT_FALSE(gmx_ana_index_has_complete_elems(&g_, INDEX_RES, top));
+
+ setGroup(group4);
+ EXPECT_FALSE(gmx_ana_index_has_complete_elems(&g_, INDEX_RES, top));
}
TEST_F(IndexBlockTest, ChecksGroupForCompleteMoleculesPositive)
topManager_.initAtoms(15);
topManager_.initUniformMolecules(3);
- t_topology *top = topManager_.topology();
+ gmx_mtop_t *top = topManager_.topology();
setGroup(group);
EXPECT_TRUE(gmx_ana_index_has_complete_elems(&g_, INDEX_MOL, top));
topManager_.initAtoms(18);
topManager_.initUniformMolecules(3);
- t_topology *top = topManager_.topology();
+ gmx_mtop_t *top = topManager_.topology();
setGroup(group1);
EXPECT_FALSE(gmx_ana_index_has_complete_elems(&g_, INDEX_MOL, top));
void PositionCalculationTest::generateCoordinates()
{
- t_topology *top = topManager_.topology();
+ t_atoms &atoms = topManager_.atoms();
t_trxframe *frame = topManager_.frame();
- for (int i = 0; i < top->atoms.nr; ++i)
+ for (int i = 0; i < atoms.nr; ++i)
{
frame->x[i][XX] = i;
- frame->x[i][YY] = top->atoms.atom[i].resind;
+ frame->x[i][YY] = atoms.atom[i].resind;
frame->x[i][ZZ] = 0.0;
if (frame->bV)
{
void
SelectionCollectionTest::setTopology()
{
- ASSERT_NO_THROW_GMX(sc_.setTopology(topManager_.mtop(), -1));
+ ASSERT_NO_THROW_GMX(sc_.setTopology(topManager_.topology(), -1));
}
void
runTest("simple.pdb", selections);
}
-TEST_F(SelectionCollectionDataTest, DISABLED_HandlesMolIndex)
+TEST_F(SelectionCollectionDataTest, HandlesMolIndex)
{
static const char * const selections[] = {
"molindex 1 4",
"molecule 2 3 5"
};
ASSERT_NO_FATAL_FAILURE(runParser(selections));
- ASSERT_NO_FATAL_FAILURE(loadTopology("simple.gro"));
+ ASSERT_NO_FATAL_FAILURE(topManager_.loadTopology("simple.gro"));
topManager_.initUniformMolecules(3);
+ ASSERT_NO_FATAL_FAILURE(setTopology());
ASSERT_NO_FATAL_FAILURE(runCompiler());
}
{
topManager_.loadTopology(filename);
- ASSERT_NO_THROW_GMX(sc_.setTopology(topManager_.mtop(), -1));
+ ASSERT_NO_THROW_GMX(sc_.setTopology(topManager_.topology(), -1));
}
#include "gromacs/fileio/trxio.h"
#include "gromacs/math/vec.h"
#include "gromacs/topology/atoms.h"
+#include "gromacs/topology/mtop_util.h"
#include "gromacs/topology/topology.h"
#include "gromacs/trajectory/trajectoryframe.h"
#include "gromacs/utility/arrayref.h"
{
TopologyManager::TopologyManager()
- : mtop_(nullptr), top_(nullptr), frame_(nullptr)
+ : mtop_(nullptr), frame_(nullptr)
{
}
done_mtop(mtop_, TRUE);
sfree(mtop_);
}
- else if (top_ != nullptr)
- {
- done_top(top_);
- }
- sfree(top_);
if (frame_ != nullptr)
{
void TopologyManager::requestFrame()
{
- GMX_RELEASE_ASSERT(top_ == NULL,
+ GMX_RELEASE_ASSERT(mtop_ == NULL,
"Frame must be requested before initializing topology");
if (frame_ == NULL)
{
void TopologyManager::initAtoms(int count)
{
- GMX_RELEASE_ASSERT(top_ == NULL, "Topology initialized more than once");
- snew(top_, 1);
- init_t_atoms(&top_->atoms, count, FALSE);
+ GMX_RELEASE_ASSERT(mtop_ == NULL, "Topology initialized more than once");
+ snew(mtop_, 1);
+ mtop_->nmoltype = 1;
+ snew(mtop_->moltype, 1);
+ init_t_atoms(&mtop_->moltype[0].atoms, count, FALSE);
+ mtop_->nmolblock = 1;
+ snew(mtop_->molblock, 1);
+ mtop_->molblock[0].type = 0;
+ mtop_->molblock[0].nmol = 1;
+ mtop_->molblock[0].natoms_mol = count;
+ mtop_->natoms = count;
+ mtop_->maxres_renum = 0;
+ gmx_mtop_finalize(mtop_);
+ GMX_RELEASE_ASSERT(mtop_->maxres_renum == 0, "maxres_renum in mtop can be modified by an env.var., that is not supported in this test");
+ t_atoms &atoms = this->atoms();
for (int i = 0; i < count; ++i)
{
- top_->atoms.atom[i].m = (i % 3 == 0 ? 2.0 : 1.0);
+ atoms.atom[i].m = (i % 3 == 0 ? 2.0 : 1.0);
}
- top_->atoms.haveMass = TRUE;
+ atoms.haveMass = TRUE;
if (frame_ != NULL)
{
frame_->natoms = count;
void TopologyManager::initUniformResidues(int residueSize)
{
- GMX_RELEASE_ASSERT(top_ != NULL, "Topology not initialized");
- int residueIndex = -1;
- for (int i = 0; i < top_->atoms.nr; ++i)
+ GMX_RELEASE_ASSERT(mtop_ != NULL, "Topology not initialized");
+ t_atoms &atoms = this->atoms();
+ int residueIndex = -1;
+ for (int i = 0; i < atoms.nr; ++i)
{
if (i % residueSize == 0)
{
++residueIndex;
}
- top_->atoms.atom[i].resind = residueIndex;
+ atoms.atom[i].resind = residueIndex;
}
}
void TopologyManager::initUniformMolecules(int moleculeSize)
{
- GMX_RELEASE_ASSERT(top_ != NULL, "Topology not initialized");
+ GMX_RELEASE_ASSERT(mtop_ != NULL, "Topology not initialized");
int index = 0;
- top_->mols.nalloc_index = (top_->atoms.nr + moleculeSize - 1) / moleculeSize + 1;
- snew(top_->mols.index, top_->mols.nalloc_index);
- top_->mols.nr = 0;
- while (index < top_->atoms.nr)
+ mtop_->mols.nalloc_index = (mtop_->natoms + moleculeSize - 1) / moleculeSize + 1;
+ snew(mtop_->mols.index, mtop_->mols.nalloc_index);
+ mtop_->mols.nr = 0;
+ while (index < mtop_->natoms)
{
- top_->mols.index[top_->mols.nr] = index;
- ++top_->mols.nr;
+ mtop_->mols.index[mtop_->mols.nr] = index;
+ ++mtop_->mols.nr;
index += moleculeSize;
}
- top_->mols.index[top_->mols.nr] = top_->atoms.nr;
+ mtop_->mols.index[mtop_->mols.nr] = mtop_->natoms;
}
void TopologyManager::initFrameIndices(const ConstArrayRef<int> &index)
struct gmx_mtop_t;
struct t_atoms;
-struct t_topology;
struct t_trxframe;
namespace gmx
void initFrameIndices(const ConstArrayRef<int> &index);
- t_topology *topology() { return top_; }
- gmx_mtop_t *mtop() { return mtop_; }
+ gmx_mtop_t *topology() { return mtop_; }
t_atoms &atoms();
t_trxframe *frame() { return frame_; }
private:
gmx_mtop_t *mtop_;
- t_topology *top_;
t_trxframe *frame_;
std::vector<char *> atomtypes_;
};
done_blocka(&(top->excls));
}
+bool gmx_mtop_has_masses(const gmx_mtop_t *mtop)
+{
+ if (mtop == nullptr)
+ {
+ return false;
+ }
+ return mtop->nmoltype == 0 || mtop->moltype[0].atoms.haveMass;
+}
+
+bool gmx_mtop_has_charges(const gmx_mtop_t *mtop)
+{
+ if (mtop == nullptr)
+ {
+ return false;
+ }
+ return mtop->nmoltype == 0 || mtop->moltype[0].atoms.haveCharge;
+}
+
static void pr_grps(FILE *fp, const char *title, const t_grps grps[], char **grpname[])
{
int i, j;
void done_mtop(gmx_mtop_t *mtop, gmx_bool bDoneSymtab);
void done_top(t_topology *top);
+bool gmx_mtop_has_masses(const gmx_mtop_t *mtop);
+bool gmx_mtop_has_charges(const gmx_mtop_t *mtop);
+
void pr_mtop(FILE *fp, int indent, const char *title, const gmx_mtop_t *mtop,
gmx_bool bShowNumbers);
void pr_top(FILE *fp, int indent, const char *title, const t_topology *top, gmx_bool bShowNumbers);
//! Returns true if a full topology file was loaded.
bool hasFullTopology() const { return bTop_; }
//! Returns the loaded topology, or NULL if not loaded.
+ const gmx_mtop_t *mtop() const { return mtop_; }
+ //! Returns the loaded topology, or NULL if not loaded.
t_topology *topology() const;
//! Returns the ePBC field from the topology.
int ePBC() const { return ePBC_; }
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2016, 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.
}
//! Helper function to initialize the grouping for a selection.
-int initSelectionGroups(Selection *sel, t_topology *top, int type)
+int initSelectionGroups(Selection *sel, const gmx_mtop_t *top, int type)
{
e_index_t indexType = INDEX_UNKNOWN;
switch (type)
PairDistance::initAnalysis(const TrajectoryAnalysisSettings &settings,
const TopologyInformation &top)
{
- refGroupCount_ = initSelectionGroups(&refSel_, top.topology(), refGroupType_);
+ refGroupCount_ = initSelectionGroups(&refSel_, top.mtop(), refGroupType_);
maxGroupCount_ = 0;
distances_.setDataSetCount(sel_.size());
for (size_t i = 0; i < sel_.size(); ++i)
{
const int selGroupCount
- = initSelectionGroups(&sel_[i], top.topology(), selGroupType_);
+ = initSelectionGroups(&sel_[i], top.mtop(), selGroupType_);
const int columnCount = refGroupCount_ * selGroupCount;
maxGroupCount_ = std::max(maxGroupCount_, columnCount);
distances_.setColumnCount(i, columnCount);
GMX_THROW(InconsistentInputError("-surf only works with -ref that consists of atoms"));
}
const e_index_t type = (surface_ == SurfaceType_Molecule ? INDEX_MOL : INDEX_RES);
- surfaceGroupCount_ = refSel_.initOriginalIdsToGroup(top.topology(), type);
+ surfaceGroupCount_ = refSel_.initOriginalIdsToGroup(top.mtop(), type);
}
if (bExclusions_)
dgsFactor_.reserve(surfaceSel_.posCount());
}
- const int resCount = surfaceSel_.initOriginalIdsToGroup(top_, INDEX_RES);
+ const int resCount = surfaceSel_.initOriginalIdsToGroup(top.mtop(), INDEX_RES);
// TODO: Not exception-safe, but nice solution would be to have a C++
// atom properties class...
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff