From: Christian Blau Date: Tue, 7 Jul 2020 18:33:21 +0000 (+0000) Subject: Move finalising routine into mtop X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=c22968c7b882512ecc45bbbb7e01de7026b841ed;p=alexxy%2Fgromacs.git Move finalising routine into mtop Change-Id: I383ceb3b48e5eddafc3f8ea96fc3b326d34602ef --- diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp index e207690390..c762178c8b 100644 --- a/src/gromacs/fileio/tpxio.cpp +++ b/src/gromacs/fileio/tpxio.cpp @@ -2837,7 +2837,7 @@ static void do_tpx_mtop(gmx::ISerializer* serializer, TpxFileHeader* tpx, gmx_mt { do_mtop(serializer, mtop, tpx->fileVersion); set_disres_npair(mtop); - gmx_mtop_finalize(mtop); + mtop->finalize(); } else { diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index 960e3e196b..7e74e01660 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -647,7 +647,7 @@ static void new_status(const char* topfile, /* Copy structures from msys to sys */ molinfo2mtop(*mi, sys); - gmx_mtop_finalize(sys); + sys->finalize(); /* Discourage using the, previously common, setup of applying constraints * to all bonds with force fields that have been parametrized with H-bond diff --git a/src/gromacs/mdlib/membed.cpp b/src/gromacs/mdlib/membed.cpp index 06266039e6..a32c758b39 100644 --- a/src/gromacs/mdlib/membed.cpp +++ b/src/gromacs/mdlib/membed.cpp @@ -1373,7 +1373,7 @@ gmx_membed_t* init_membed(FILE* fplog, // Re-establish the invariants of the derived values within // mtop. - gmx_mtop_finalize(mtop); + mtop->finalize(); if (ftp2bSet(efTOP, nfile, fnm)) { diff --git a/src/gromacs/mdlib/tests/updategroupscog.cpp b/src/gromacs/mdlib/tests/updategroupscog.cpp index ab21b65fd4..b77e92b359 100644 --- a/src/gromacs/mdlib/tests/updategroupscog.cpp +++ b/src/gromacs/mdlib/tests/updategroupscog.cpp @@ -1,7 +1,7 @@ /* * This file is part of the GROMACS molecular simulation package. * - * Copyright (c) 2018,2019, by the GROMACS development team, led by + * Copyright (c) 2018,2019,2020, 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. @@ -98,7 +98,7 @@ TEST(UpdateGroupsCog, ComputesCogs) mtop.molblock[0].type = 0; mtop.molblock[0].nmol = numMolecules; mtop.natoms = numAtoms; - gmx_mtop_finalize(&mtop); + mtop.finalize(); // Set up the SETTLE parameters. const real dOH = 0.1; diff --git a/src/gromacs/pbcutil/tests/com.cpp b/src/gromacs/pbcutil/tests/com.cpp index 02e5b87a82..84249c46b1 100644 --- a/src/gromacs/pbcutil/tests/com.cpp +++ b/src/gromacs/pbcutil/tests/com.cpp @@ -153,7 +153,7 @@ COMInPlaceTest::COMInPlaceTest() : molblock.nmol = 2; molblock.type = 0; testTopology_.natoms = moltype.atoms.nr * molblock.nmol; - gmx_mtop_finalize(&testTopology_); + testTopology_.finalize(); FloatingPointTolerance tolerance( FloatingPointTolerance(1.0e-6, 1.0e-6, 1.0e-8, 1.0e-12, 10000, 100, false)); checker_.setDefaultTolerance(tolerance); diff --git a/src/gromacs/selection/tests/toputils.cpp b/src/gromacs/selection/tests/toputils.cpp index 77885f7913..8b6388db44 100644 --- a/src/gromacs/selection/tests/toputils.cpp +++ b/src/gromacs/selection/tests/toputils.cpp @@ -153,9 +153,8 @@ void TopologyManager::initAtoms(int count) mtop_->molblock[0].type = 0; mtop_->molblock[0].nmol = 1; mtop_->natoms = count; - mtop_->maxres_renum = 0; - gmx_mtop_finalize(mtop_.get()); - GMX_RELEASE_ASSERT(mtop_->maxres_renum == 0, + mtop_->finalize(); + GMX_RELEASE_ASSERT(mtop_->maxResiduesPerMoleculeToTriggerRenumber() == 0, "maxres_renum in mtop can be modified by an env.var., that is not supported " "in this test"); t_atoms& atoms = this->atoms(); @@ -254,9 +253,8 @@ void TopologyManager::finalizeTopology() { GMX_RELEASE_ASSERT(mtop_ != nullptr, "Topology not initialized"); - mtop_->maxres_renum = 0; mtop_->haveMoleculeIndices = true; - gmx_mtop_finalize(mtop_.get()); + mtop_->finalize(); } void TopologyManager::initUniformResidues(int residueSize) @@ -291,7 +289,7 @@ void TopologyManager::initUniformMolecules(int moleculeSize) "The residues should break at molecule boundaries"); atoms.nres = nres; mtop_->haveMoleculeIndices = true; - gmx_mtop_finalize(mtop_.get()); + mtop_->finalize(); } void TopologyManager::initFrameIndices(const ArrayRef& index) diff --git a/src/gromacs/topology/mtop_lookup.h b/src/gromacs/topology/mtop_lookup.h index f6491bbc9c..0c61cc0265 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,2018,2019, by the GROMACS development team, led by + * Copyright (c) 2016,2017,2018,2019,2020, 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. @@ -219,7 +219,7 @@ static inline void mtopGetAtomAndResidueName(const gmx_mtop_t* mtop, } if (residueNumber != nullptr) { - if (atoms.nres > mtop->maxres_renum) + if (atoms.nres > mtop->maxResiduesPerMoleculeToTriggerRenumber()) { *residueNumber = atoms.resinfo[atoms.atom[atomIndexInMolecule].resind].nr; } diff --git a/src/gromacs/topology/mtop_util.cpp b/src/gromacs/topology/mtop_util.cpp index 0c2006def8..fb4a6ea7d6 100644 --- a/src/gromacs/topology/mtop_util.cpp +++ b/src/gromacs/topology/mtop_util.cpp @@ -57,94 +57,6 @@ #include "gromacs/utility/real.h" #include "gromacs/utility/smalloc.h" -static int gmx_mtop_maxresnr(const gmx_mtop_t* mtop, int maxres_renum) -{ - int maxresnr = 0; - - for (const gmx_moltype_t& moltype : mtop->moltype) - { - const t_atoms& atoms = moltype.atoms; - if (atoms.nres > maxres_renum) - { - for (int r = 0; r < atoms.nres; r++) - { - if (atoms.resinfo[r].nr > maxresnr) - { - maxresnr = atoms.resinfo[r].nr; - } - } - } - } - - return maxresnr; -} - -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 (size_t mb = 0; mb < mtop->molblock.size(); mb++) - { - const gmx_molblock_t& molb = mtop->molblock[mb]; - MoleculeBlockIndices& indices = mtop->moleculeBlockIndices[mb]; - const int numResPerMol = mtop->moltype[molb.type].atoms.nres; - - indices.numAtomsPerMolecule = mtop->moltype[molb.type].atoms.nr; - indices.globalAtomStart = atomIndex; - indices.globalResidueStart = residueIndex; - atomIndex += molb.nmol * indices.numAtomsPerMolecule; - residueIndex += molb.nmol * numResPerMol; - indices.globalAtomEnd = atomIndex; - indices.residueNumberStart = residueNumberStart; - if (numResPerMol <= mtop->maxres_renum) - { - residueNumberStart += molb.nmol * numResPerMol; - } - indices.moleculeIndexStart = moleculeIndexStart; - moleculeIndexStart += molb.nmol; - } -} - -void gmx_mtop_finalize(gmx_mtop_t* mtop) -{ - char* env; - - if (mtop->molblock.size() == 1 && mtop->molblock[0].nmol == 1) - { - /* We have a single molecule only, no renumbering needed. - * This case also covers an mtop converted from pdb/gro/... input, - * so we retain the original residue numbering. - */ - mtop->maxres_renum = 0; - } - else - { - /* We only renumber single residue molecules. Their intra-molecular - * residue numbering is anyhow irrelevant. - */ - mtop->maxres_renum = 1; - } - - env = getenv("GMX_MAXRESRENUM"); - if (env != nullptr) - { - sscanf(env, "%d", &mtop->maxres_renum); - } - if (mtop->maxres_renum == -1) - { - /* -1 signals renumber residues in all molecules */ - mtop->maxres_renum = INT_MAX; - } - - mtop->maxresnr = gmx_mtop_maxresnr(mtop, mtop->maxres_renum); - - buildMolblockIndices(mtop); -} - void gmx_mtop_count_atomtypes(const gmx_mtop_t* mtop, int state, int typecount[]) { for (int i = 0; i < mtop->ffparams.atnr; ++i) @@ -195,7 +107,7 @@ AtomIterator::AtomIterator(const gmx_mtop_t& mtop, int globalAtomNumber) : mblock_(0), atoms_(&mtop.moltype[mtop.molblock[0].type].atoms), currentMolecule_(0), - highestResidueNumber_(mtop.maxresnr), + highestResidueNumber_(mtop.maxResNumberNotRenumbered()), localAtomNumber_(0), globalAtomNumber_(globalAtomNumber) { @@ -210,7 +122,7 @@ AtomIterator& AtomIterator::operator++() if (localAtomNumber_ >= atoms_->nr) { - if (atoms_->nres <= mtop_->maxresnr) + if (atoms_->nres <= mtop_->maxResiduesPerMoleculeToTriggerRenumber()) { /* Single residue molecule, increase the count with one */ highestResidueNumber_ += atoms_->nres; @@ -260,7 +172,7 @@ const char* AtomProxy::residueName() const int AtomProxy::residueNumber() const { int residueIndexInMolecule = it_->atoms_->atom[it_->localAtomNumber_].resind; - if (it_->atoms_->nres <= it_->mtop_->maxres_renum) + if (it_->atoms_->nres <= it_->mtop_->maxResiduesPerMoleculeToTriggerRenumber()) { return it_->highestResidueNumber_ + 1 + residueIndexInMolecule; } @@ -580,10 +492,11 @@ t_atoms gmx_mtop_global_atoms(const gmx_mtop_t* mtop) init_t_atoms(&atoms, 0, FALSE); - int maxresnr = mtop->maxresnr; + int maxresnr = mtop->maxResNumberNotRenumbered(); for (const gmx_molblock_t& molb : mtop->molblock) { - atomcat(&atoms, &mtop->moltype[molb.type].atoms, molb.nmol, mtop->maxres_renum, &maxresnr); + atomcat(&atoms, &mtop->moltype[molb.type].atoms, molb.nmol, + mtop->maxResiduesPerMoleculeToTriggerRenumber(), &maxresnr); } return atoms; @@ -1135,7 +1048,7 @@ void convertAtomsToMtop(t_symtab* symtab, char** name, t_atoms* atoms, gmx_mtop_ mtop->haveMoleculeIndices = false; - gmx_mtop_finalize(mtop); + mtop->finalize(); } bool haveFepPerturbedNBInteractions(const gmx_mtop_t* mtop) diff --git a/src/gromacs/topology/mtop_util.h b/src/gromacs/topology/mtop_util.h index e027bfef7c..9a027bac7d 100644 --- a/src/gromacs/topology/mtop_util.h +++ b/src/gromacs/topology/mtop_util.h @@ -58,12 +58,6 @@ struct t_symtab; // and should be replaced by versions taking const gmx_mtop & when // their callers are refactored similarly. -/* Should be called after generating or reading mtop, - * to set some compute intesive variables to avoid - * N^2 operations later on. - */ -void gmx_mtop_finalize(gmx_mtop_t* mtop); - /* Counts the number of atoms of each type. State should be 0 for * state A and 1 for state B types. typecount should have at * least mtop->ffparams.atnr elements. diff --git a/src/gromacs/topology/tests/mtop.cpp b/src/gromacs/topology/tests/mtop.cpp index 517ab10ba0..08713b39eb 100644 --- a/src/gromacs/topology/tests/mtop.cpp +++ b/src/gromacs/topology/tests/mtop.cpp @@ -68,7 +68,7 @@ void createBasicTop(gmx_mtop_t* mtop) mtop->molblock[0].type = 0; mtop->molblock[0].nmol = 3; mtop->natoms = moltype.atoms.nr * mtop->molblock[0].nmol; - gmx_mtop_finalize(mtop); + mtop->finalize(); } /*! \brief @@ -99,7 +99,7 @@ std::vector> createTwoResidueTopology(gmx_mtop_t* mtop) mtop->molblock[0].type = 0; mtop->molblock[0].nmol = 1; mtop->natoms = moltype.atoms.nr * mtop->molblock[0].nmol; - gmx_mtop_finalize(mtop); + mtop->finalize(); std::vector> residueRange; residueRange.emplace_back(0, residueOneSize); residueRange.emplace_back(residueOneSize, residueOneSize + residueTwoSize); diff --git a/src/gromacs/topology/topology.cpp b/src/gromacs/topology/topology.cpp index f24c788d4f..6c9bfc1d60 100644 --- a/src/gromacs/topology/topology.cpp +++ b/src/gromacs/topology/topology.cpp @@ -86,6 +86,28 @@ gmx_moltype_t::~gmx_moltype_t() done_atom(&atoms); } +static int gmx_mtop_maxresnr(const gmx::ArrayRef moltypes, int maxres_renum) +{ + int maxresnr = 0; + + for (const gmx_moltype_t& moltype : moltypes) + { + const t_atoms& atoms = moltype.atoms; + if (atoms.nres > maxres_renum) + { + for (int r = 0; r < atoms.nres; r++) + { + if (atoms.resinfo[r].nr > maxresnr) + { + maxresnr = atoms.resinfo[r].nr; + } + } + } + } + + return maxresnr; +} + gmx_mtop_t::gmx_mtop_t() { init_atomtypes(&atomtypes); @@ -101,6 +123,70 @@ gmx_mtop_t::~gmx_mtop_t() done_atomtypes(&atomtypes); } +void gmx_mtop_t::finalize() +{ + if (molblock.size() == 1 && molblock[0].nmol == 1) + { + /* We have a single molecule only, no renumbering needed. + * This case also covers an mtop converted from pdb/gro/... input, + * so we retain the original residue numbering. + */ + maxResiduesPerMoleculeToTriggerRenumber_ = 0; + } + else + { + /* We only renumber single residue molecules. Their intra-molecular + * residue numbering is anyhow irrelevant. + */ + maxResiduesPerMoleculeToTriggerRenumber_ = 1; + } + + const char* env = getenv("GMX_MAXRESRENUM"); + if (env != nullptr) + { + sscanf(env, "%d", &maxResiduesPerMoleculeToTriggerRenumber_); + } + if (maxResiduesPerMoleculeToTriggerRenumber_ == -1) + { + /* -1 signals renumber residues in all molecules */ + maxResiduesPerMoleculeToTriggerRenumber_ = std::numeric_limits::max(); + } + + maxResNumberNotRenumbered_ = gmx_mtop_maxresnr(moltype, maxResiduesPerMoleculeToTriggerRenumber_); + + buildMolblockIndices(); +} + +void gmx_mtop_t::buildMolblockIndices() +{ + moleculeBlockIndices.resize(molblock.size()); + + int atomIndex = 0; + int residueIndex = 0; + int residueNumberStart = maxResNumberNotRenumbered_ + 1; + int moleculeIndexStart = 0; + for (size_t mb = 0; mb < molblock.size(); mb++) + { + const gmx_molblock_t& molb = molblock[mb]; + MoleculeBlockIndices& indices = moleculeBlockIndices[mb]; + const int numResPerMol = moltype[molb.type].atoms.nres; + + indices.numAtomsPerMolecule = moltype[molb.type].atoms.nr; + indices.globalAtomStart = atomIndex; + indices.globalResidueStart = residueIndex; + atomIndex += molb.nmol * indices.numAtomsPerMolecule; + residueIndex += molb.nmol * numResPerMol; + indices.globalAtomEnd = atomIndex; + indices.residueNumberStart = residueNumberStart; + if (numResPerMol <= maxResiduesPerMoleculeToTriggerRenumber_) + { + residueNumberStart += molb.nmol * numResPerMol; + } + indices.moleculeIndexStart = moleculeIndexStart; + moleculeIndexStart += molb.nmol; + } +} + void done_top(t_topology* top) { done_idef(&top->idef); @@ -593,8 +679,9 @@ void compareMtop(FILE* fp, const gmx_mtop_t& mtop1, const gmx_mtop_t& mtop2, rea fprintf(fp, "comparing mtop topology\n"); cmp_str(fp, "Name", -1, *mtop1.name, *mtop2.name); cmp_int(fp, "natoms", -1, mtop1.natoms, mtop2.natoms); - cmp_int(fp, "maxres_renum", -1, mtop1.maxres_renum, mtop2.maxres_renum); - cmp_int(fp, "maxresnr", -1, mtop1.maxresnr, mtop2.maxresnr); + cmp_int(fp, "maxres_renum", -1, mtop1.maxResiduesPerMoleculeToTriggerRenumber(), + mtop2.maxResiduesPerMoleculeToTriggerRenumber()); + cmp_int(fp, "maxresnr", -1, mtop1.maxResNumberNotRenumbered(), mtop2.maxResNumberNotRenumbered()); cmp_bool(fp, "bIntermolecularInteractions", -1, mtop1.bIntermolecularInteractions, mtop2.bIntermolecularInteractions); cmp_bool(fp, "haveMoleculeIndices", -1, mtop1.haveMoleculeIndices, mtop2.haveMoleculeIndices); diff --git a/src/gromacs/topology/topology.h b/src/gromacs/topology/topology.h index a3fe344f77..6ac65d0b39 100644 --- a/src/gromacs/topology/topology.h +++ b/src/gromacs/topology/topology.h @@ -177,10 +177,6 @@ struct gmx_mtop_t //NOLINT(clang-analyzer-optin.performance.Padding) std::unique_ptr intermolecular_ilist = nullptr; //! Number of global atoms. int natoms = 0; - //! Parameter for residue numbering. - int maxres_renum = 0; - //! The maximum residue number in moltype - int maxresnr = -1; //! Atomtype properties t_atomtypes atomtypes; //! Groups of atoms for different purposes @@ -194,9 +190,33 @@ struct gmx_mtop_t //NOLINT(clang-analyzer-optin.performance.Padding) */ std::vector intermolecularExclusionGroup; - /* Derived data below */ + //! Maximum number of residues in molecule to trigger renumbering of residues + int maxResiduesPerMoleculeToTriggerRenumber() const + { + return maxResiduesPerMoleculeToTriggerRenumber_; + } + //! Maximum residue number that is not renumbered. + int maxResNumberNotRenumbered() const { return maxResNumberNotRenumbered_; } + /*! \brief Finalize this data structure. + * + * Should be called after generating or reading mtop, to set some compute + * intesive variables to avoid N^2 operations later on. + * + * \todo Move into a builder class, once available. + */ + void finalize(); + + /* Derived data below */ //! Indices for each molblock entry for fast lookup of atom properties std::vector moleculeBlockIndices; + +private: + //! Build the molblock indices + void buildMolblockIndices(); + //! Maximum number of residues in molecule to trigger renumbering of residues + int maxResiduesPerMoleculeToTriggerRenumber_ = 0; + //! The maximum residue number in moltype that is not renumbered + int maxResNumberNotRenumbered_ = -1; }; /*! \brief