From 8dd3c9ae88004054b3b112c23f747d27a19d8d29 Mon Sep 17 00:00:00 2001 From: Berk Hess Date: Wed, 14 Feb 2018 15:36:49 +0100 Subject: [PATCH] Remove mols from gmx_mtop_t The molecule atom begin/end information in gmx_mtop_t was derived information. Now the mols struct is removed and the information is derived when necessary. Introduced BlockRanges, a C++ verion of t_block. This is preparation for converison of gmx_mtop_t to C++. Change-Id: I1717a101de82354a7f9cf14fd7f43c7227e554ef --- src/gromacs/domdec/domdec.cpp | 2 +- src/gromacs/fileio/tpxio.cpp | 35 +-------------- src/gromacs/gmxana/gmx_clustsize.cpp | 17 ++++---- src/gromacs/imd/imd.cpp | 8 ++-- src/gromacs/mdlib/broadcaststructs.cpp | 4 +- src/gromacs/mdlib/forcerec.cpp | 3 +- src/gromacs/mdlib/tests/settle.cpp | 1 - src/gromacs/selection/indexutil.cpp | 47 +++++++++++--------- src/gromacs/selection/indexutil.h | 4 +- src/gromacs/selection/sm_simple.cpp | 17 +++----- src/gromacs/selection/tests/indexutil.cpp | 7 ++- src/gromacs/selection/tests/toputils.cpp | 27 ++++++------ src/gromacs/tools/check.cpp | 2 +- src/gromacs/topology/block.h | 45 ++++++++++++++++---- src/gromacs/topology/mtop_lookup.h | 25 ++++++++++- src/gromacs/topology/mtop_util.cpp | 71 ++++++++++++++++++++++++++++++- src/gromacs/topology/mtop_util.h | 18 +++++++- src/gromacs/topology/topology.cpp | 5 +-- src/gromacs/topology/topology.h | 5 ++- src/programs/mdrun/membed.cpp | 36 ++++++---------- 20 files changed, 240 insertions(+), 139 deletions(-) diff --git a/src/gromacs/domdec/domdec.cpp b/src/gromacs/domdec/domdec.cpp index eb5305dafd..801fe558eb 100644 --- a/src/gromacs/domdec/domdec.cpp +++ b/src/gromacs/domdec/domdec.cpp @@ -6384,7 +6384,7 @@ static void set_dd_limits_and_grid(FILE *fplog, t_commrec *cr, gmx_domdec_t *dd, comm->bCGs = (ncg_mtop(mtop) < mtop->natoms); - comm->bInterCGBondeds = ((ncg_mtop(mtop) > mtop->mols.nr) || + comm->bInterCGBondeds = ((ncg_mtop(mtop) > gmx_mtop_num_molecules(*mtop)) || mtop->bIntermolecularInteractions); if (comm->bInterCGBondeds) { diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp index d60a5c97ed..7a841ec380 100644 --- a/src/gromacs/fileio/tpxio.cpp +++ b/src/gromacs/fileio/tpxio.cpp @@ -2445,35 +2445,6 @@ static void do_molblock(t_fileio *fio, gmx_molblock_t *molb, gmx_bool bRead) } -static t_block mtop_mols(gmx_mtop_t *mtop) -{ - int mb, m, a, mol; - t_block mols; - - mols.nr = 0; - for (mb = 0; mb < mtop->nmolblock; mb++) - { - mols.nr += mtop->molblock[mb].nmol; - } - mols.nalloc_index = mols.nr + 1; - snew(mols.index, mols.nalloc_index); - - a = 0; - m = 0; - mols.index[m] = a; - for (mb = 0; mb < mtop->nmolblock; mb++) - { - for (mol = 0; mol < mtop->molblock[mb].nmol; mol++) - { - a += mtop->molblock[mb].natoms_mol; - m++; - mols.index[m] = a; - } - } - - return mols; -} - static void set_disres_npair(gmx_mtop_t *mtop) { t_iparams *ip; @@ -2575,11 +2546,7 @@ static void do_mtop(t_fileio *fio, gmx_mtop_t *mtop, gmx_bool bRead, do_groups(fio, &mtop->groups, bRead, &(mtop->symtab)); - if (bRead) - { - done_block(&mtop->mols); - mtop->mols = mtop_mols(mtop); - } + mtop->haveMoleculeIndices = true; if (bRead) { diff --git a/src/gromacs/gmxana/gmx_clustsize.cpp b/src/gromacs/gmxana/gmx_clustsize.cpp index 738843ac99..5ea9d9fa79 100644 --- a/src/gromacs/gmxana/gmx_clustsize.cpp +++ b/src/gromacs/gmxana/gmx_clustsize.cpp @@ -54,6 +54,7 @@ #include "gromacs/pbcutil/pbc.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/trajectory/trajectoryframe.h" #include "gromacs/utility/arraysize.h" @@ -83,7 +84,6 @@ static void clust_size(const char *ndx, const char *trx, const char *xpm, t_tpxheader tpxh; gmx_mtop_t *mtop = nullptr; int ePBC = -1; - t_block *mols = nullptr; int ii, jj; real temp, tfac; /* Cluster size distribution (matrix) */ @@ -132,6 +132,7 @@ static void clust_size(const char *ndx, const char *trx, const char *xpm, tfac = ndf/(3.0*natoms); } + gmx::BlockRanges mols; if (bMol) { if (ndx) @@ -140,10 +141,10 @@ static void clust_size(const char *ndx, const char *trx, const char *xpm, ndx); } GMX_RELEASE_ASSERT(mtop != nullptr, "Trying to access mtop->mols from NULL mtop pointer"); - mols = &(mtop->mols); + mols = gmx_mtop_molecules(*mtop); /* Make dummy index */ - nindex = mols->nr; + nindex = mols.numBlocks(); snew(index, nindex); for (i = 0; (i < nindex); i++) { @@ -209,11 +210,11 @@ static void clust_size(const char *ndx, const char *trx, const char *xpm, /* Compute distance */ if (bMol) { - GMX_RELEASE_ASSERT(mols != nullptr, "Cannot access index[] from NULL mols pointer"); + GMX_RELEASE_ASSERT(mols.numBlocks() > 0, "Cannot access index[] from empty mols"); bSame = FALSE; - for (ii = mols->index[ai]; !bSame && (ii < mols->index[ai+1]); ii++) + for (ii = mols.index[ai]; !bSame && (ii < mols.index[ai+1]); ii++) { - for (jj = mols->index[aj]; !bSame && (jj < mols->index[aj+1]); jj++) + for (jj = mols.index[aj]; !bSame && (jj < mols.index[aj+1]); jj++) { if (bPBC) { @@ -364,8 +365,8 @@ static void clust_size(const char *ndx, const char *trx, const char *xpm, { if (bMol) { - GMX_RELEASE_ASSERT(mols != nullptr, "Cannot access index[] from NULL mols pointer"); - for (j = mols->index[i]; (j < mols->index[i+1]); j++) + GMX_RELEASE_ASSERT(mols.numBlocks() > 0, "Cannot access index[] from empty mols"); + for (j = mols.index[i]; (j < mols.index[i+1]); j++) { fprintf(fp, "%d\n", j+1); } diff --git a/src/gromacs/imd/imd.cpp b/src/gromacs/imd/imd.cpp index 02666d3c80..bc242a4f34 100644 --- a/src/gromacs/imd/imd.cpp +++ b/src/gromacs/imd/imd.cpp @@ -1048,11 +1048,10 @@ static void init_imd_prepare_mols_in_imdgroup(t_gmx_IMD_setup *IMDsetup, gmx_mto { int i, ii; int gstart, gend, count; - t_block gmols, lmols; + t_block lmols; int nat; int *ind; - gmols = top_global->mols; nat = IMDsetup->nat; ind = IMDsetup->ind; @@ -1067,10 +1066,11 @@ static void init_imd_prepare_mols_in_imdgroup(t_gmx_IMD_setup *IMDsetup, gmx_mto } } - snew(lmols.index, gmols.nr+1); + gmx::BlockRanges gmols = gmx_mtop_molecules(*top_global); + snew(lmols.index, gmols.numBlocks() + 1); lmols.index[0] = 0; - for (i = 0; i < gmols.nr; i++) + for (i = 0; i < gmols.numBlocks(); i++) { gstart = gmols.index[i]; gend = gmols.index[i+1]; diff --git a/src/gromacs/mdlib/broadcaststructs.cpp b/src/gromacs/mdlib/broadcaststructs.cpp index 6a570aacc5..e47995f19d 100644 --- a/src/gromacs/mdlib/broadcaststructs.cpp +++ b/src/gromacs/mdlib/broadcaststructs.cpp @@ -825,8 +825,10 @@ void bcast_ir_mtop(const t_commrec *cr, t_inputrec *inputrec, gmx_mtop_t *mtop) bc_atomtypes(cr, &mtop->atomtypes); - bc_block(cr, &mtop->mols); 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; } void init_parallel(t_commrec *cr, t_inputrec *inputrec, diff --git a/src/gromacs/mdlib/forcerec.cpp b/src/gromacs/mdlib/forcerec.cpp index b77cf9a602..d2618016fd 100644 --- a/src/gromacs/mdlib/forcerec.cpp +++ b/src/gromacs/mdlib/forcerec.cpp @@ -2358,7 +2358,8 @@ void init_forcerec(FILE *fp, */ cgs = &mtop->moltype[mtop->molblock[mtop->nmolblock-1].type].cgs; fr->n_tpi = cgs->index[cgs->nr] - cgs->index[cgs->nr-1]; - if (fr->n_tpi != mtop->mols.index[mtop->mols.nr] - mtop->mols.index[mtop->mols.nr-1]) + gmx::BlockRanges molecules = gmx_mtop_molecules(*mtop); + if (fr->n_tpi != molecules.index[molecules.numBlocks()] - molecules.index[molecules.numBlocks() - 1]) { gmx_fatal(FARGS, "The molecule to insert can not consist of multiple charge groups.\nMake it a single charge group."); } diff --git a/src/gromacs/mdlib/tests/settle.cpp b/src/gromacs/mdlib/tests/settle.cpp index f605b0d581..f717fc0207 100644 --- a/src/gromacs/mdlib/tests/settle.cpp +++ b/src/gromacs/mdlib/tests/settle.cpp @@ -201,7 +201,6 @@ TEST_P(SettleTest, SatisfiesConstraints) gmx_mtop_t *mtop; snew(mtop, 1); const unique_cptr mtopGuard(mtop); - mtop->mols.nr = 1; mtop->nmoltype = 1; snew(mtop->moltype, mtop->nmoltype); const unique_cptr moltypeGuard(mtop->moltype); diff --git a/src/gromacs/selection/indexutil.cpp b/src/gromacs/selection/indexutil.cpp index d264b47bb6..b2c1b5b353 100644 --- a/src/gromacs/selection/indexutil.cpp +++ b/src/gromacs/selection/indexutil.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016,2017, by the GROMACS development team, led by + * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016,2017,2018, 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. @@ -796,7 +796,8 @@ gmx_ana_index_merge(gmx_ana_index_t *dest, * \ingroup module_selection */ static bool -next_group_index(int atomIndex, const gmx_mtop_t *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) @@ -812,16 +813,11 @@ next_group_index(int atomIndex, const gmx_mtop_t *top, e_index_t type, int *id) break; } case INDEX_MOL: - if (*id >= 0 && top->mols.index[*id] > atomIndex) - { - *id = 0; - } - while (*id < top->mols.nr && atomIndex >= top->mols.index[*id+1]) - { - ++*id; - } - GMX_ASSERT(*id < top->mols.nr, "Molecules do not span all the atoms"); + { + int molb = 0; + *id = mtopGetMoleculeIndex(top, atomIndex, &molb); break; + } case INDEX_UNKNOWN: case INDEX_ALL: *id = 0; @@ -868,7 +864,7 @@ gmx_ana_index_make_block(t_blocka *t, const gmx_mtop_t *top, gmx_ana_index_t *g, // into exceptions. GMX_RELEASE_ASSERT(top != nullptr || (type != INDEX_RES && type != INDEX_MOL), "Topology must be provided for residue or molecule blocks"); - GMX_RELEASE_ASSERT(!(type == INDEX_MOL && top->mols.nr == 0), + GMX_RELEASE_ASSERT(type != INDEX_MOL || top->haveMoleculeIndices, "Molecule information must be present for molecule blocks"); /* bComplete only does something for INDEX_RES or INDEX_MOL, so turn it @@ -956,12 +952,20 @@ gmx_ana_index_make_block(t_blocka *t, const gmx_mtop_t *top, gmx_ana_index_t *g, break; } case INDEX_MOL: - for (int j = top->mols.index[id]; j < top->mols.index[id+1]; ++j) + { + int molb = 0; + while (molb + 1 < top->nmolblock && id >= top->molblock[molb].moleculeIndexStart) { - t->a[t->nra++] = j; + ++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) + { + t->a[t->nra++] = atomStart + j; } break; - + } default: /* Should not be reached */ GMX_RELEASE_ASSERT(false, "Unreachable code was reached"); break; @@ -995,7 +999,7 @@ gmx_ana_index_make_block(t_blocka *t, const gmx_mtop_t *top, gmx_ana_index_t *g, * The atoms in \p g are assumed to be sorted. */ bool -gmx_ana_index_has_full_blocks(const gmx_ana_index_t *g, const t_block *b) +gmx_ana_index_has_full_blocks(const gmx_ana_index_t *g, const gmx::BlockRanges *b) { int i, j, bi; @@ -1004,12 +1008,12 @@ gmx_ana_index_has_full_blocks(const gmx_ana_index_t *g, const t_block *b) while (i < g->isize) { /* Find the block that begins with the first unmatched atom */ - while (bi < b->nr && b->index[bi] != g->index[i]) + while (bi < b->numBlocks() && b->index[bi] != g->index[i]) { ++bi; } /* If not found, or if too large, return */ - if (bi == b->nr || i + b->index[bi+1] - b->index[bi] > g->isize) + if (bi == b->numBlocks() || i + b->index[bi+1] - b->index[bi] > g->isize) { return false; } @@ -1154,7 +1158,10 @@ gmx_ana_index_has_complete_elems(gmx_ana_index_t *g, e_index_t type, } case INDEX_MOL: - return gmx_ana_index_has_full_blocks(g, &top->mols); + { + gmx::BlockRanges molecules = gmx_mtop_molecules(*top); + return gmx_ana_index_has_full_blocks(g, &molecules); + } } return true; } @@ -1257,8 +1264,6 @@ gmx_ana_indexmap_init_orgid_group(gmx_ana_indexmap_t *m, const gmx_mtop_t *top, "to use the mapping"); GMX_RELEASE_ASSERT(top != nullptr || (type != INDEX_RES && type != INDEX_MOL), "Topology must be provided for residue or molecule blocks"); - GMX_RELEASE_ASSERT(!(type == INDEX_MOL && top->mols.nr == 0), - "Molecule information must be present for molecule blocks"); // Check that all atoms in each block belong to the same group. // This is a separate loop for better error handling (no state is modified // if there is an error. diff --git a/src/gromacs/selection/indexutil.h b/src/gromacs/selection/indexutil.h index 83de646fef..099ad91ed1 100644 --- a/src/gromacs/selection/indexutil.h +++ b/src/gromacs/selection/indexutil.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016, by the GROMACS development team, led by + * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016,2018, 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. @@ -331,7 +331,7 @@ 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(const gmx_ana_index_t *g, const t_block *b); +gmx_ana_index_has_full_blocks(const gmx_ana_index_t *g, const gmx::BlockRanges *b); /** Checks whether a group consists of full blocks. */ bool gmx_ana_index_has_full_ablocks(gmx_ana_index_t *g, t_blocka *b); diff --git a/src/gromacs/selection/sm_simple.cpp b/src/gromacs/selection/sm_simple.cpp index 259e395fe6..619a19f2c2 100644 --- a/src/gromacs/selection/sm_simple.cpp +++ b/src/gromacs/selection/sm_simple.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016,2017, by the GROMACS development team, led by + * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016,2017,2018, 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. @@ -601,7 +601,7 @@ check_molecules(const gmx_mtop_t *top, int /* npar */, gmx_ana_selparam_t * /* p { bool bOk; - bOk = (top != nullptr && top->mols.nr > 0); + bOk = (top != nullptr && top->haveMoleculeIndices); if (!bOk) { GMX_THROW(gmx::InconsistentInputError("Molecule information not available in topology")); @@ -618,16 +618,11 @@ static void evaluate_molindex(const gmx::SelMethodEvalContext &context, gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void * /* data */) { - int i, j; - - out->nr = g->isize; - for (i = j = 0; i < g->isize; ++i) + out->nr = g->isize; + int molb = 0; + for (int i = 0; i < g->isize; ++i) { - while (context.top->mols.index[j + 1] <= g->index[i]) - { - ++j; - } - out->u.i[i] = j + 1; + out->u.i[i] = mtopGetMoleculeIndex(context.top, g->index[i], &molb) + 1; } } diff --git a/src/gromacs/selection/tests/indexutil.cpp b/src/gromacs/selection/tests/indexutil.cpp index a2d43aca86..c91519e343 100644 --- a/src/gromacs/selection/tests/indexutil.cpp +++ b/src/gromacs/selection/tests/indexutil.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017,2018, 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. @@ -177,6 +177,7 @@ TEST_F(IndexBlockTest, CreatesMoleculeBlock) { const int group[] = { 3, 4, 7, 8, 13 }; topManager_.initAtoms(18); + topManager_.initUniformResidues(1); topManager_.initUniformMolecules(3); setGroup(group); gmx_ana_index_make_block(&blocka_, topManager_.topology(), &g_, @@ -199,6 +200,7 @@ TEST_F(IndexBlockTest, CreatesMoleculeBlockWithCompletion) { const int group[] = { 3, 4, 7, 8, 13 }; topManager_.initAtoms(18); + topManager_.initUniformResidues(1); topManager_.initUniformMolecules(3); setGroup(group); gmx_ana_index_make_block(&blocka_, topManager_.topology(), &g_, @@ -334,6 +336,7 @@ TEST_F(IndexBlockTest, ChecksGroupForCompleteMoleculesPositive) const int group[] = { 0, 1, 2, 6, 7, 8, 12, 13, 14 }; topManager_.initAtoms(15); + topManager_.initUniformResidues(1); topManager_.initUniformMolecules(3); gmx_mtop_t *top = topManager_.topology(); @@ -348,6 +351,7 @@ TEST_F(IndexBlockTest, ChecksGroupForCompleteMoleculesNegative) const int group3[] = { 4, 5, 6, 7, 8, 12, 13, 14 }; topManager_.initAtoms(18); + topManager_.initUniformResidues(1); topManager_.initUniformMolecules(3); gmx_mtop_t *top = topManager_.topology(); @@ -531,6 +535,7 @@ TEST_F(IndexMapTest, InitializesMoleculeBlock) { const int maxGroup[] = { 3, 4, 7, 8, 13 }; topManager_.initAtoms(18); + topManager_.initUniformResidues(1); topManager_.initUniformMolecules(3); testInit(maxGroup, INDEX_MOL); } diff --git a/src/gromacs/selection/tests/toputils.cpp b/src/gromacs/selection/tests/toputils.cpp index 79eecdf31a..234f466309 100644 --- a/src/gromacs/selection/tests/toputils.cpp +++ b/src/gromacs/selection/tests/toputils.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017,2018, 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. @@ -226,22 +226,25 @@ void TopologyManager::initUniformResidues(int residueSize) } atoms.atom[i].resind = residueIndex; } + atoms.nres = residueIndex; } void TopologyManager::initUniformMolecules(int moleculeSize) { GMX_RELEASE_ASSERT(mtop_ != nullptr, "Topology not initialized"); - int index = 0; - mtop_->mols.nalloc_index = (mtop_->natoms + moleculeSize - 1) / moleculeSize + 1; - srenew(mtop_->mols.index, mtop_->mols.nalloc_index); - mtop_->mols.nr = 0; - while (index < mtop_->natoms) - { - mtop_->mols.index[mtop_->mols.nr] = index; - ++mtop_->mols.nr; - index += moleculeSize; - } - mtop_->mols.index[mtop_->mols.nr] = mtop_->natoms; + GMX_RELEASE_ASSERT(mtop_->nmolblock == 1, "initUniformMolecules only implemented for a single molblock"); + gmx_molblock_t &molblock = mtop_->molblock[0]; + t_atoms &atoms = mtop_->moltype[molblock.type].atoms; + GMX_RELEASE_ASSERT(atoms.nr % moleculeSize == 0, + "The number of atoms should be a multiple of moleculeSize"); + molblock.nmol = atoms.nr/moleculeSize; + atoms.nr = moleculeSize; + const int nres = atoms.atom[atoms.nr].resind; + GMX_RELEASE_ASSERT(atoms.atom[atoms.nr-1].resind != nres, + "The residues should break at molecule boundaries"); + atoms.nres = nres; + molblock.natoms_mol = moleculeSize; + mtop_->haveMoleculeIndices = true; } void TopologyManager::initFrameIndices(const ArrayRef &index) diff --git a/src/gromacs/tools/check.cpp b/src/gromacs/tools/check.cpp index d7429f6a2c..e9ab81b00f 100644 --- a/src/gromacs/tools/check.cpp +++ b/src/gromacs/tools/check.cpp @@ -204,7 +204,7 @@ static void tpx2system(FILE *fp, const gmx_mtop_t *mtop) } } fprintf(fp, "A system of %d molecules (%d atoms) was simulated.\n", - mtop->mols.nr, mtop->natoms-nvsite); + gmx_mtop_num_molecules(*mtop), mtop->natoms-nvsite); if (nvsite) { fprintf(fp, "Virtual sites were used in some of the molecules.\n"); diff --git a/src/gromacs/topology/block.h b/src/gromacs/topology/block.h index d4bbe40c13..0131f185ba 100644 --- a/src/gromacs/topology/block.h +++ b/src/gromacs/topology/block.h @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2010,2014,2015, by the GROMACS development team, led by + * Copyright (c) 2010,2014,2015,2018, 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. @@ -39,20 +39,49 @@ #include +#ifdef __cplusplus +#include +#endif + #include "gromacs/utility/basedefinitions.h" +#ifdef __cplusplus +#include "gromacs/utility/gmxassert.h" +#endif #ifdef __cplusplus extern "C" { #endif -/* the block structure points into an array (usually of ints). - It is a list of starting indices for objects of consecutive ids, such - as molecules. - For example, if this block denotes molecules, then the first molecule - ranges from index[0] to index[1]-1 in the atom list. +#ifdef __cplusplus +namespace gmx +{ + +/*! \brief Division of a range of indices into consecutive blocks + * + * A range of consecutive indices 0 to index[numBlocks()] is divided + * into numBlocks() consecutive blocks of consecutive indices. + * Block b contains indices i for which index[b] <= i < index[b+1]. + */ +struct BlockRanges +{ + /*! \brief Returns the number of blocks + * + * This should only be called on a valid struct. Validy is asserted + * (only) in debug mode. + */ + int numBlocks() const + { + GMX_ASSERT(index.size() > 0, "numBlocks() should only be called on a valid BlockRanges struct"); + return index.size() - 1; + } + + std::vector index; /**< The list of block begin/end indices */ +}; + +} // nsamespace gmx +#endif // __cplusplus - This makes the mapping from atoms to molecules O(Nmolecules) instead - of O(Natoms) in size. */ +/* Deprecated, C-style version of BlockRanges */ typedef struct t_block { int nr; /* The number of blocks */ diff --git a/src/gromacs/topology/mtop_lookup.h b/src/gromacs/topology/mtop_lookup.h index 1d892ab179..70853e0e70 100644 --- a/src/gromacs/topology/mtop_lookup.h +++ b/src/gromacs/topology/mtop_lookup.h @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2016,2017, by the GROMACS development team, led by + * Copyright (c) 2016,2017,2018, 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. @@ -110,6 +110,29 @@ mtopGetMolblockIndex(const gmx_mtop_t *mtop, } } +/*! \brief Returns the global molecule index of a global atom index + * + * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms. + * The input value of moleculeBlock should be in range. Use 0 as starting value. + * For subsequent calls to this function, e.g. in a loop, pass in the previously + * returned value for best performance. Atoms in a group tend to be in the same + * molecule(block), so this minimizes the search time. + * + * \param[in] mtop The molecule topology + * \param[in] globalAtomIndex The global atom index to look up + * \param[in,out] moleculeBlock The molecule block index in \p mtop + */ +static inline int +mtopGetMoleculeIndex(const gmx_mtop_t *mtop, + int globalAtomIndex, + int *moleculeBlock) +{ + int localMoleculeIndex; + mtopGetMolblockIndex(mtop, globalAtomIndex, moleculeBlock, &localMoleculeIndex, nullptr); + + return mtop->molblock[*moleculeBlock].moleculeIndexStart + localMoleculeIndex; +} + /*! \brief Returns the atom data for an atom based on global atom index * * The atom index has to be in range: 0 <= \p globalAtomIndex < \p mtop->natoms. diff --git a/src/gromacs/topology/mtop_util.cpp b/src/gromacs/topology/mtop_util.cpp index 52b1305bd3..96dfa454b2 100644 --- a/src/gromacs/topology/mtop_util.cpp +++ b/src/gromacs/topology/mtop_util.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2008,2009,2010,2012,2013,2014,2015,2016,2017, by the GROMACS development team, led by + * Copyright (c) 2008,2009,2010,2012,2013,2014,2015,2016,2017,2018, 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. @@ -49,6 +49,7 @@ #include "gromacs/topology/ifunc.h" #include "gromacs/topology/topology.h" #include "gromacs/topology/topsort.h" +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/real.h" #include "gromacs/utility/smalloc.h" @@ -83,6 +84,7 @@ static void finalizeMolblocks(gmx_mtop_t *mtop) int atomIndex = 0; int residueIndex = 0; int residueNumberStart = mtop->maxresnr + 1; + int moleculeIndexStart = 0; for (int mb = 0; mb < mtop->nmolblock; mb++) { gmx_molblock_t *molb = &mtop->molblock[mb]; @@ -97,6 +99,8 @@ static void finalizeMolblocks(gmx_mtop_t *mtop) { residueNumberStart += molb->nmol*numResPerMol; } + molb->moleculeIndexStart = moleculeIndexStart; + moleculeIndexStart += molb->nmol; } } @@ -180,6 +184,16 @@ int ncg_mtop(const gmx_mtop_t *mtop) return ncg; } +int gmx_mtop_num_molecules(const gmx_mtop_t &mtop) +{ + int numMolecules = 0; + for (int mb = 0; mb < mtop.nmolblock; mb++) + { + numMolecules += mtop.molblock[mb].nmol; + } + return numMolecules; +} + int gmx_mtop_nres(const gmx_mtop_t *mtop) { int nres = 0; @@ -968,6 +982,57 @@ gmx_mtop_generate_local_top(const gmx_mtop_t *mtop, return top; } +/*! \brief Fills an array with molecule begin/end atom indices + * + * \param[in] mtop The global topology + * \param[out] index Array of size nr. of molecules + 1 to be filled with molecule begin/end indices + */ +static void fillMoleculeIndices(const gmx_mtop_t &mtop, + gmx::ArrayRef index) +{ + int globalAtomIndex = 0; + int globalMolIndex = 0; + index[globalMolIndex] = globalAtomIndex; + for (int mb = 0; mb < mtop.nmolblock; mb++) + { + const gmx_molblock_t &molb = mtop.molblock[mb]; + for (int mol = 0; mol < molb.nmol; mol++) + { + globalAtomIndex += molb.natoms_mol; + globalMolIndex += 1; + index[globalMolIndex] = globalAtomIndex; + } + } +} + +gmx::BlockRanges gmx_mtop_molecules(const gmx_mtop_t &mtop) +{ + gmx::BlockRanges mols; + + mols.index.resize(gmx_mtop_num_molecules(mtop) + 1); + + fillMoleculeIndices(mtop, mols.index); + + return mols; +} + +/*! \brief Creates and returns a deprecated t_block struct with molecule indices + * + * \param[in] mtop The global topology + */ +static t_block gmx_mtop_molecules_t_block(const gmx_mtop_t &mtop) +{ + t_block mols; + + mols.nr = gmx_mtop_num_molecules(mtop); + mols.nalloc_index = mols.nr + 1; + snew(mols.index, mols.nalloc_index); + + fillMoleculeIndices(mtop, gmx::arrayRefFromArray(mols.index, mols.nr + 1)); + + return mols; +} + t_topology gmx_mtop_t_to_t_topology(gmx_mtop_t *mtop, bool freeMTop) { int mt, mb; @@ -983,7 +1048,7 @@ t_topology gmx_mtop_t_to_t_topology(gmx_mtop_t *mtop, bool freeMTop) top.cgs = ltop.cgs; top.excls = ltop.excls; top.atoms = gmx_mtop_global_atoms(mtop); - top.mols = mtop->mols; + top.mols = gmx_mtop_molecules_t_block(*mtop); top.bIntermolecularInteractions = mtop->bIntermolecularInteractions; top.symtab = mtop->symtab; @@ -1053,5 +1118,7 @@ void convertAtomsToMtop(t_symtab *symtab, mtop->natoms = atoms->nr; + mtop->haveMoleculeIndices = false; + gmx_mtop_finalize(mtop); } diff --git a/src/gromacs/topology/mtop_util.h b/src/gromacs/topology/mtop_util.h index 366cd13e9e..918bb3dcfd 100644 --- a/src/gromacs/topology/mtop_util.h +++ b/src/gromacs/topology/mtop_util.h @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2012,2013,2014,2015,2016, by the GROMACS development team, led by + * Copyright (c) 2012,2013,2014,2015,2016,2018, 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. @@ -70,6 +70,12 @@ gmx_mtop_count_atomtypes(const gmx_mtop_t *mtop, int state, int typecount[]); int ncg_mtop(const gmx_mtop_t *mtop); +/*!\brief Returns the total number of molecules in mtop + * + * \param[in] mtop The global topology + */ +int gmx_mtop_num_molecules(const gmx_mtop_t &mtop); + /* Returns the total number of residues in mtop. */ int gmx_mtop_nres(const gmx_mtop_t *mtop); @@ -210,6 +216,16 @@ gmx_localtop_t * gmx_mtop_generate_local_top(const gmx_mtop_t *mtop, bool freeEnergyInteractionsAtEnd); +/*!\brief Creates and returns a struct with begin/end atom indices of all molecules + * + * \param[in] mtop The global topology + * \returns A struct of type BlockRanges with numBlocks() equal to the number + * of molecules and atom indices such that molecule m contains atoms a with: + * index[m] <= a < index[m+1]. + */ +gmx::BlockRanges gmx_mtop_molecules(const gmx_mtop_t &mtop); + + /* Converts a gmx_mtop_t struct to t_topology. * * If freeMTop == true, memory related to mtop will be freed so that done_top() diff --git a/src/gromacs/topology/topology.cpp b/src/gromacs/topology/topology.cpp index ac9ec0d9bb..a9b6b312eb 100644 --- a/src/gromacs/topology/topology.cpp +++ b/src/gromacs/topology/topology.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017,2018, 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. @@ -80,7 +80,6 @@ void init_mtop(gmx_mtop_t *mtop) mtop->maxres_renum = 0; mtop->maxresnr = -1; init_groups(&mtop->groups); - init_block(&mtop->mols); open_symtab(&mtop->symtab); } @@ -168,7 +167,6 @@ void done_mtop(gmx_mtop_t *mtop) sfree(mtop->molblock); done_atomtypes(&mtop->atomtypes); done_gmx_groups_t(&mtop->groups); - done_block(&mtop->mols); } void done_top(t_topology *top) @@ -208,6 +206,7 @@ void done_top_mtop(t_topology *top, gmx_mtop_t *mtop) done_atom(&top->atoms); done_block(&top->cgs); done_blocka(&top->excls); + done_block(&top->mols); done_symtab(&top->symtab); open_symtab(&mtop->symtab); } diff --git a/src/gromacs/topology/topology.h b/src/gromacs/topology/topology.h index 89c40bfc08..c9a226bd3e 100644 --- a/src/gromacs/topology/topology.h +++ b/src/gromacs/topology/topology.h @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2011,2014,2015,2016, by the GROMACS development team, led by + * Copyright (c) 2011,2014,2015,2016,2018, 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. @@ -79,6 +79,7 @@ typedef struct gmx_molblock_t 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 */ int residueNumberStart; /**< Residue numbers start from this value if the number of residues per molecule is <= maxres_renum */ + int moleculeIndexStart; /**< Global molecule indexing starts from this value */ } gmx_molblock_t; typedef struct gmx_groups_t @@ -115,9 +116,9 @@ typedef struct gmx_mtop_t int maxres_renum; /* Parameter for residue numbering */ int maxresnr; /* The maximum residue number in moltype */ t_atomtypes atomtypes; /* Atomtype properties */ - t_block mols; /* The molecules */ gmx_groups_t groups; t_symtab symtab; /* The symbol table */ + bool haveMoleculeIndices; /* Tells whether we have valid molecule indices */ } gmx_mtop_t; /* The mdrun node-local topology struct, completely written out */ diff --git a/src/programs/mdrun/membed.cpp b/src/programs/mdrun/membed.cpp index 90db647cb0..e6f09b7380 100644 --- a/src/programs/mdrun/membed.cpp +++ b/src/programs/mdrun/membed.cpp @@ -206,8 +206,8 @@ static void check_types(t_block *ins_at, t_block *rest_at, gmx_mtop_t *mtop) } } - sfree(ins_mtype->index); - sfree(rest_mtype->index); + done_block(ins_mtype); + done_block(rest_mtype); sfree(ins_mtype); sfree(rest_mtype); } @@ -546,8 +546,8 @@ static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc int *order; r_min_rad = probe_rad*probe_rad; - snew(rm_p->mol, mtop->mols.nr); - snew(rm_p->block, mtop->mols.nr); + gmx::BlockRanges molecules = gmx_mtop_molecules(*mtop); + snew(rm_p->block, molecules.numBlocks()); nrm = nupper = 0; nlower = low_up_rm; for (i = 0; i < ins_at->nr; i++) @@ -580,7 +580,7 @@ static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc { if (mol_id == mem_p->mol_id[l]) { - for (k = mtop->mols.index[mol_id]; k < mtop->mols.index[mol_id+1]; k++) + for (k = molecules.index[mol_id]; k < molecules.index[mol_id+1]; k++) { z_lip += r[k][ZZ]; } @@ -607,7 +607,7 @@ static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc snew(order, mem_p->nmol); for (i = 0; i < mem_p->nmol; i++) { - at = mtop->mols.index[mem_p->mol_id[i]]; + at = molecules.index[mem_p->mol_id[i]]; pbc_dx(pbc, r[at], pos_ins->geom_cent[0], dr); if (pos_ins->pieces > 1) { @@ -650,7 +650,7 @@ static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc if (bRM) { z_lip = 0; - for (k = mtop->mols.index[mol_id]; k < mtop->mols.index[mol_id+1]; k++) + for (k = molecules.index[mol_id]; k < molecules.index[mol_id+1]; k++) { z_lip += r[k][ZZ]; } @@ -694,17 +694,20 @@ static void rm_group(gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state { int i, j, k, n, rm, mol_id, at, block; rvec *x_tmp, *v_tmp; - int *list, *new_mols; + int *list; unsigned char *new_egrp[egcNR]; gmx_bool bRM; int RMmolblock; + /* Construct the molecule range information */ + gmx::BlockRanges molecules = gmx_mtop_molecules(*mtop); + snew(list, state->natoms); n = 0; for (i = 0; i < rm_p->nr; i++) { mol_id = rm_p->mol[i]; - at = mtop->mols.index[mol_id]; + at = molecules.index[mol_id]; block = rm_p->block[i]; mtop->molblock[block].nmol--; for (j = 0; j < mtop->molblock[block].natoms_mol; j++) @@ -712,23 +715,8 @@ static void rm_group(gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state list[n] = at+j; n++; } - mtop->mols.index[mol_id] = -1; } - mtop->mols.nr -= rm_p->nr; - mtop->mols.nalloc_index -= rm_p->nr; - snew(new_mols, mtop->mols.nr); - for (i = 0; i < mtop->mols.nr+rm_p->nr; i++) - { - j = 0; - if (mtop->mols.index[i] != -1) - { - new_mols[j] = mtop->mols.index[i]; - j++; - } - } - sfree(mtop->mols.index); - mtop->mols.index = new_mols; mtop->natoms -= n; state->natoms -= n; snew(x_tmp, state->natoms); -- 2.11.4.GIT