Changed posres pointers in molblock to std vector.
Also moved the derived helper index data into a separate struct.
That is merged into this change to keep MPI broadcasting simple.
natoms_mol will be moved into the new struct in a child change.
Change-Id: I131fd14ab53cec98800c7e0731d3a92df3f93422
/* Get the position restraint coordinates from the molblock */
a_molb = mol*molb->natoms_mol + a_mol;
- if (a_molb >= molb->nposres_xA)
+ if (a_molb >= static_cast<int>(molb->posres_xA.size()))
{
gmx_incons("Not enough position restraint coordinates");
}
ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
- if (molb->nposres_xB > 0)
+ if (!molb->posres_xB.empty())
{
ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
/* Get the position restriant coordinats from the molblock */
a_molb = mol*molb->natoms_mol + a_mol;
- if (a_molb >= molb->nposres_xA)
+ if (a_molb >= static_cast<int>(molb->posres_xA.size()))
{
gmx_incons("Not enough position restraint coordinates");
}
gmx_fio_do_int(fio, molb->nmol);
gmx_fio_do_int(fio, molb->natoms_mol);
/* Position restraint coordinates */
- gmx_fio_do_int(fio, molb->nposres_xA);
- if (molb->nposres_xA > 0)
+ int numPosres_xA = molb->posres_xA.size();
+ gmx_fio_do_int(fio, numPosres_xA);
+ if (numPosres_xA > 0)
{
if (bRead)
{
- snew(molb->posres_xA, molb->nposres_xA);
+ molb->posres_xA.resize(numPosres_xA);
}
- gmx_fio_ndo_rvec(fio, molb->posres_xA, molb->nposres_xA);
+ gmx_fio_ndo_rvec(fio, as_rvec_array(molb->posres_xA.data()), numPosres_xA);
}
- gmx_fio_do_int(fio, molb->nposres_xB);
- if (molb->nposres_xB > 0)
+ int numPosres_xB = molb->posres_xB.size();
+ gmx_fio_do_int(fio, numPosres_xB);
+ if (numPosres_xB > 0)
{
if (bRead)
{
- snew(molb->posres_xB, molb->nposres_xB);
+ molb->posres_xB.resize(numPosres_xB);
}
- gmx_fio_ndo_rvec(fio, molb->posres_xB, molb->nposres_xB);
+ gmx_fio_ndo_rvec(fio, as_rvec_array(molb->posres_xB.data()), numPosres_xB);
}
}
sys->molblock.push_back(molb);
gmx_molblock_t &molbs = sys->molblock.back();
molbs.natoms_mol = molinfo[molbs.type].atoms.nr;
- molbs.nposres_xA = 0;
- molbs.nposres_xB = 0;
sys->natoms += molbs.nmol*molbs.natoms_mol;
}
}
warninp_t wi)
{
gmx_bool *hadAtom;
- rvec *x, *v, *xp;
+ rvec *x, *v;
dvec sum;
double totmass;
t_topology *top;
}
if (!bTopB)
{
- molb.nposres_xA = nat_molb;
- snew(molb.posres_xA, molb.nposres_xA);
+ molb.posres_xA.resize(nat_molb);
for (i = 0; i < nat_molb; i++)
{
copy_rvec(x[a+i], molb.posres_xA[i]);
}
else
{
- molb.nposres_xB = nat_molb;
- snew(molb.posres_xB, molb.nposres_xB);
+ molb.posres_xB.resize(nat_molb);
for (i = 0; i < nat_molb; i++)
{
copy_rvec(x[a+i], molb.posres_xB[i]);
for (gmx_molblock_t &molb : mtop->molblock)
{
nat_molb = molb.nmol*mtop->moltype[molb.type].atoms.nr;
- if (molb.nposres_xA > 0 || molb.nposres_xB > 0)
+ if (!molb.posres_xA.empty() || !molb.posres_xB.empty())
{
- xp = (!bTopB ? molb.posres_xA : molb.posres_xB);
+ std::vector<gmx::RVec> &xp = (!bTopB ? molb.posres_xA : molb.posres_xB);
for (i = 0; i < nat_molb; i++)
{
for (j = 0; j < npbcdim; j++)
bc_blocka(cr, &moltype->excls);
}
-static void bc_molblock(const t_commrec *cr, gmx_molblock_t *molb)
+static void bc_vector_of_rvec(const t_commrec *cr, std::vector<gmx::RVec> *vec)
{
- block_bc(cr, *molb);
- if (molb->nposres_xA > 0)
+ int numElements = vec->size();
+ block_bc(cr, numElements);
+ if (!MASTER(cr))
{
- snew_bc(cr, molb->posres_xA, molb->nposres_xA);
- nblock_bc(cr, molb->nposres_xA*DIM, molb->posres_xA[0]);
+ vec->resize(numElements);
}
- if (molb->nposres_xB > 0)
+ if (numElements > 0)
{
- snew_bc(cr, molb->posres_xB, molb->nposres_xB);
- nblock_bc(cr, molb->nposres_xB*DIM, molb->posres_xB[0]);
+ nblock_bc(cr, numElements, as_rvec_array(vec->data()));
}
+}
+
+static void bc_molblock(const t_commrec *cr, gmx_molblock_t *molb)
+{
+ block_bc(cr, molb->type);
+ block_bc(cr, molb->nmol);
+ bc_vector_of_rvec(cr, &molb->posres_xA);
+ bc_vector_of_rvec(cr, &molb->posres_xB);
+ block_bc(cr, molb->natoms_mol);
if (debug)
{
fprintf(debug, "after bc_molblock\n");
bc_groups(cr, &mtop->symtab, mtop->natoms, &mtop->groups);
GMX_RELEASE_ASSERT(!MASTER(cr) || mtop->haveMoleculeIndices, "mtop should have valid molecule indices");
- mtop->haveMoleculeIndices = true;
+ if (!MASTER(cr))
+ {
+ mtop->haveMoleculeIndices = true;
+
+ gmx_mtop_finalize(mtop);
+ }
}
void init_parallel(t_commrec *cr, t_inputrec *inputrec,
gmx::ArrayRef<gmx::RVec> x)
{
GMX_ASSERT(x.size() >= static_cast<size_t>(mtop.natoms), "x should contain the whole system");
+ GMX_ASSERT(!mtop.moleculeBlockIndices.empty(), "molblock indices are needed in constructVsitesGlobal");
- for (const gmx_molblock_t &molb : mtop.molblock)
+ for (size_t mb = 0; mb < mtop.molblock.size(); mb++)
{
- const gmx_moltype_t &molt = mtop.moltype[molb.type];
+ const gmx_molblock_t &molb = mtop.molblock[mb];
+ const gmx_moltype_t &molt = mtop.moltype[molb.type];
if (vsiteIlistNrCount(molt.ilist) > 0)
{
- int atomOffset = molb.globalAtomStart;
+ int atomOffset = mtop.moleculeBlockIndices[mb].globalAtomStart;
for (int mol = 0; mol < molb.nmol; mol++)
{
construct_vsites(nullptr, as_rvec_array(x.data()) + atomOffset,
{
--first_atom;
}
- int first_mol_atom = top->molblock[molb].globalAtomStart;
+ int first_mol_atom = top->moleculeBlockIndices[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;
case INDEX_MOL:
{
size_t molb = 0;
- while (molb + 1 < top->molblock.size() && id >= top->molblock[molb].moleculeIndexStart)
+ while (molb + 1 < top->molblock.size() && id >= top->moleculeBlockIndices[molb].moleculeIndexStart)
{
++molb;
}
- const gmx_molblock_t &molblock = top->molblock[molb];
- const int atomStart = molblock.globalAtomStart + (id - molblock.moleculeIndexStart)*molblock.natoms_mol;
- for (int j = 0; j < molblock.natoms_mol; ++j)
+ const int numAtomsInMol = top->molblock[molb].natoms_mol;
+ const MoleculeBlockIndices &blockIndices = top->moleculeBlockIndices[molb];
+ const int atomStart = blockIndices.globalAtomStart + (id - blockIndices.moleculeIndexStart)*numAtomsInMol;
+ for (int j = 0; j < numAtomsInMol; ++j)
{
t->a[t->nra++] = atomStart + j;
}
mtop->molblock[0].type = 0;
mtop->molblock[0].nmol = 1;
mtop->molblock[0].natoms_mol = top.atoms.nr;
- mtop->molblock[0].nposres_xA = 0;
- mtop->molblock[0].nposres_xB = 0;
mtop->natoms = top.atoms.nr;
}
GMX_ASSERT(globalAtomIndex >= 0 && globalAtomIndex < mtop->natoms, "The atom index to look up should be within range");
GMX_ASSERT(moleculeBlock != nullptr, "molBlock can not be NULL");
GMX_ASSERT(*moleculeBlock >= 0 && *moleculeBlock < static_cast<int>(mtop->molblock.size()), "The starting molecule block index for the search should be within range");
+ GMX_ASSERT(!mtop->moleculeBlockIndices.empty(), "The molecule block indices should be present when calling mtopGetMoleculeBlockIndex");
/* Search the molecue block index using bisection */
int molBlock0 = -1;
int globalAtomStart;
while (TRUE)
{
- globalAtomStart = mtop->molblock[*moleculeBlock].globalAtomStart;
+ globalAtomStart = mtop->moleculeBlockIndices[*moleculeBlock].globalAtomStart;
if (globalAtomIndex < globalAtomStart)
{
molBlock1 = *moleculeBlock;
}
- else if (globalAtomIndex >= mtop->molblock[*moleculeBlock].globalAtomEnd)
+ else if (globalAtomIndex >= mtop->moleculeBlockIndices[*moleculeBlock].globalAtomEnd)
{
molBlock0 = *moleculeBlock;
}
int localMoleculeIndex;
mtopGetMolblockIndex(mtop, globalAtomIndex, moleculeBlock, &localMoleculeIndex, nullptr);
- return mtop->molblock[*moleculeBlock].moleculeIndexStart + localMoleculeIndex;
+ return mtop->moleculeBlockIndices[*moleculeBlock].moleculeIndexStart + localMoleculeIndex;
}
/*! \brief Returns the atom data for an atom based on global atom index
mtopGetMolblockIndex(mtop, globalAtomIndex, moleculeBlock,
&moleculeIndex, &atomIndexInMolecule);
- const gmx_molblock_t &molb = mtop->molblock[*moleculeBlock];
- const t_atoms &atoms = mtop->moltype[molb.type].atoms;
+ const gmx_molblock_t &molb = mtop->molblock[*moleculeBlock];
+ const t_atoms &atoms = mtop->moltype[molb.type].atoms;
+ const MoleculeBlockIndices &indices = mtop->moleculeBlockIndices[*moleculeBlock];
if (atomName != nullptr)
{
*atomName = *(atoms.atomname[atomIndexInMolecule]);
else
{
/* Single residue molecule, keep counting */
- *residueNumber = molb.residueNumberStart + moleculeIndex*atoms.nres + atoms.atom[atomIndexInMolecule].resind;
+ *residueNumber = indices.residueNumberStart + moleculeIndex*atoms.nres + atoms.atom[atomIndexInMolecule].resind;
}
}
if (residueName != nullptr)
}
if (globalResidueIndex != nullptr)
{
- *globalResidueIndex = molb.globalResidueStart + moleculeIndex*atoms.nres + atoms.atom[atomIndexInMolecule].resind;
+ *globalResidueIndex = indices.globalResidueStart + moleculeIndex*atoms.nres + atoms.atom[atomIndexInMolecule].resind;
}
}
return maxresnr;
}
-static void finalizeMolblocks(gmx_mtop_t *mtop)
+static void buildMolblockIndices(gmx_mtop_t *mtop)
{
+ mtop->moleculeBlockIndices.resize(mtop->molblock.size());
+
int atomIndex = 0;
int residueIndex = 0;
int residueNumberStart = mtop->maxresnr + 1;
int moleculeIndexStart = 0;
- for (gmx_molblock_t &molb : mtop->molblock)
+ for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
{
- const int numResPerMol = mtop->moltype[molb.type].atoms.nres;
- molb.globalAtomStart = atomIndex;
- molb.globalResidueStart = residueIndex;
+ const gmx_molblock_t &molb = mtop->molblock[mb];
+ MoleculeBlockIndices &indices = mtop->moleculeBlockIndices[mb];
+ const int numResPerMol = mtop->moltype[molb.type].atoms.nres;
+
+ indices.globalAtomStart = atomIndex;
+ indices.globalResidueStart = residueIndex;
atomIndex += molb.nmol*molb.natoms_mol;
residueIndex += molb.nmol*numResPerMol;
- molb.globalAtomEnd = atomIndex;
- molb.residueNumberStart = residueNumberStart;
+ indices.globalAtomEnd = atomIndex;
+ indices.residueNumberStart = residueNumberStart;
if (numResPerMol <= mtop->maxres_renum)
{
residueNumberStart += molb.nmol*numResPerMol;
}
- molb.moleculeIndexStart = moleculeIndexStart;
+ indices.moleculeIndexStart = moleculeIndexStart;
moleculeIndexStart += molb.nmol;
}
}
mtop->maxresnr = gmx_mtop_maxresnr(mtop, mtop->maxres_renum);
- finalizeMolblocks(mtop);
+ buildMolblockIndices(mtop);
}
void gmx_mtop_count_atomtypes(const gmx_mtop_t *mtop, int state, int typecount[])
/* Copy the force constants */
*ip = idef->iparams[il->iatoms[i*2]];
a_molb = il->iatoms[i*2+1] - a_offset;
- if (molb->nposres_xA == 0)
+ if (molb->posres_xA.empty())
{
gmx_incons("Position restraint coordinates are missing");
}
ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
- if (molb->nposres_xB > 0)
+ if (!molb->posres_xB.empty())
{
ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
/* Copy the force constants */
*ip = idef->iparams[il->iatoms[i*2]];
a_molb = il->iatoms[i*2+1] - a_offset;
- if (molb->nposres_xA == 0)
+ if (molb->posres_xA.empty())
{
gmx_incons("Position restraint coordinates are missing");
}
#include "gromacs/topology/ifunc.h"
#include "gromacs/topology/symtab.h"
#include "gromacs/utility/compare.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/strconvert.h"
#include "gromacs/utility/txtdump.h"
}
}
-gmx_molblock_t::gmx_molblock_t()
-{
- type = -1;
- nmol = 0;
- nposres_xA = 0;
- posres_xA = nullptr;
- nposres_xB = 0;
- posres_xB = nullptr;
-
- natoms_mol = 0;
- globalAtomStart = -1;
- globalAtomEnd = -1;
- globalResidueStart = -1;
- residueNumberStart = -1;
- moleculeIndexStart = -1;
-}
-
-gmx_molblock_t::~gmx_molblock_t()
-{
- if (nposres_xA > 0)
- {
- nposres_xA = 0;
- sfree(posres_xA);
- }
- if (nposres_xB > 0)
- {
- nposres_xB = 0;
- sfree(posres_xB);
- }
-}
-
void done_gmx_groups_t(gmx_groups_t *g)
{
int i;
"moltype", molb->type, *(molt[molb->type].name));
pr_int(fp, indent, "#molecules", molb->nmol);
pr_int(fp, indent, "#atoms_mol", molb->natoms_mol);
- pr_int(fp, indent, "#posres_xA", molb->nposres_xA);
- if (molb->nposres_xA > 0)
+ pr_int(fp, indent, "#posres_xA", molb->posres_xA.size());
+ if (!molb->posres_xA.empty())
{
- pr_rvecs(fp, indent, "posres_xA", molb->posres_xA, molb->nposres_xA);
+ pr_rvecs(fp, indent, "posres_xA", as_rvec_array(molb->posres_xA.data()), molb->posres_xA.size());
}
- pr_int(fp, indent, "#posres_xB", molb->nposres_xB);
- if (molb->nposres_xB > 0)
+ pr_int(fp, indent, "#posres_xB", molb->posres_xB.size());
+ if (!molb->posres_xB.empty())
{
- pr_rvecs(fp, indent, "posres_xB", molb->posres_xB, molb->nposres_xB);
+ pr_rvecs(fp, indent, "posres_xB", as_rvec_array(molb->posres_xB.data()), molb->posres_xB.size());
}
}
/*! \brief Block of molecules of the same type, used in gmx_mtop_t */
struct gmx_molblock_t
{
- /*! \brief Constructor */
- gmx_molblock_t();
-
- /*! \brief Destructor */
- ~gmx_molblock_t();
-
- /*! \brief Default copy assignment operator.
- *
- * NOTE: Does not free the old pointers.
- */
- gmx_molblock_t &operator=(const gmx_molblock_t &) = default;
-
- /*! \brief Default copy constructor */
- gmx_molblock_t(const gmx_molblock_t &) = default;
+ int type = -1; /**< The molecule type index in mtop.moltype */
+ int nmol = 0; /**< The number of molecules in this block */
+ std::vector<gmx::RVec> posres_xA; /**< Position restraint coordinates for top A */
+ std::vector<gmx::RVec> posres_xB; /**< Position restraint coordinates for top B */
- int type; /**< The molecule type index in mtop.moltype */
- int nmol; /**< The number of molecules in this block */
- int nposres_xA; /**< The number of posres coords for top A */
- rvec *posres_xA; /**< Position restraint coordinates for top A */
- int nposres_xB; /**< The number of posres coords for top B */
- rvec *posres_xB; /**< Position restraint coordinates for top B */
+ /* Convenience information, derived from other gmx_mtop_t contents */
+ int natoms_mol = 0; /**< The number of atoms in one molecule */
+};
- /* Convenience information, derived from other gmx_mtop_t contents */
- int natoms_mol; /**< The number of atoms in one molecule */
+/*! \brief Indices for a gmx_molblock_t, derived from other gmx_mtop_t contents */
+struct MoleculeBlockIndices
+{
int globalAtomStart; /**< Global atom index of the first atom in the block */
int globalAtomEnd; /**< Global atom index + 1 of the last atom in the block */
int globalResidueStart; /**< Global residue index of the first residue in the block */
#define ggrpnr(groups, egc, i) ((groups)->grpnr[egc] ? (groups)->grpnr[egc][i] : 0)
/* The global, complete system topology struct, based on molecule types.
- This structure should contain no data that is O(natoms) in memory. */
+ * This structure should contain no data that is O(natoms) in memory.
+ *
+ * TODO: Find a solution for ensuring that the derived data is in sync
+ * with the primary data, possibly by converting to a class.
+ */
struct gmx_mtop_t
{
/* Constructor */
int maxres_renum; /* Parameter for residue numbering */
int maxresnr; /* The maximum residue number in moltype */
t_atomtypes atomtypes; /* Atomtype properties */
- gmx_groups_t groups;
+ gmx_groups_t groups; /* Groups of atoms for different purposes */
t_symtab symtab; /* The symbol table */
bool haveMoleculeIndices; /* Tells whether we have valid molecule indices */
+
+ /* Derived data */
+ std::vector<MoleculeBlockIndices> moleculeBlockIndices; /* Indices for each molblock entry for fast lookup of atom properties */
};
/* The mdrun node-local topology struct, completely written out */
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff