Move finalising routine into mtop
authorChristian Blau <cblau.mail@gmail.com>
Tue, 7 Jul 2020 18:33:21 +0000 (18:33 +0000)
committerArtem Zhmurov <zhmurov@gmail.com>
Tue, 7 Jul 2020 18:33:21 +0000 (18:33 +0000)
Change-Id: I383ceb3b48e5eddafc3f8ea96fc3b326d34602ef

12 files changed:
src/gromacs/fileio/tpxio.cpp
src/gromacs/gmxpreprocess/grompp.cpp
src/gromacs/mdlib/membed.cpp
src/gromacs/mdlib/tests/updategroupscog.cpp
src/gromacs/pbcutil/tests/com.cpp
src/gromacs/selection/tests/toputils.cpp
src/gromacs/topology/mtop_lookup.h
src/gromacs/topology/mtop_util.cpp
src/gromacs/topology/mtop_util.h
src/gromacs/topology/tests/mtop.cpp
src/gromacs/topology/topology.cpp
src/gromacs/topology/topology.h

index e20769039028d42235c0c149972c50252aef6ff9..c762178c8b0a63d210a10d8d453fe289afe7e0ed 100644 (file)
@@ -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
         {
index 960e3e196b4275a5bc7b3170b11bb47c135ec83b..7e74e01660098e7e4e5ac0415d77e3931470a9cd 100644 (file)
@@ -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
index 06266039e6f9a7a1aafdf681f9d9b465e0e6d548..a32c758b39007ef8f5859777cb48878a95bf0a99 100644 (file)
@@ -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))
         {
index ab21b65fd4281dffd2dfee88c65cf13d4c7a6cea..b77e92b359ebfd52b89b8c3c0119b172933ca4c9 100644 (file)
@@ -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;
index 02e5b87a826d6f792fb1d397e915b071b7db0ff5..84249c46b1b2dafd7e59281f071ef3c6bdad4e44 100644 (file)
@@ -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);
index 77885f791336dc78de1de1ba24560bcdf2529af7..8b6388db44bfe893d199f9e5437902f12b53b9e3 100644 (file)
@@ -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<const int>& index)
index f6491bbc9c052cc35e716d43734ee51a36311aa4..0c61cc0265f6621558a4d8ed36407f41ec9ba597 100644 (file)
@@ -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;
         }
index 0c2006def870cacd50ec5c18ce0f75b3826db334..fb4a6ea7d615dbaf76fa20424e701b254a482d80 100644 (file)
 #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)
index e027bfef7c138a2c2056227c0885297c1d39fa59..9a027bac7dd99e970fc751c5e2fd74b46a261f6f 100644 (file)
@@ -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.
index 517ab10ba040a0cb9421016464fd6334f675705e..08713b39eb85787ef0d3f523d46169fc64464d48 100644 (file)
@@ -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<gmx::Range<int>> 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<gmx::Range<int>> residueRange;
     residueRange.emplace_back(0, residueOneSize);
     residueRange.emplace_back(residueOneSize, residueOneSize + residueTwoSize);
index f24c788d4f7d0110a827757aba65e17703005bee..6c9bfc1d60d2e8209947590076dd281aa770963d 100644 (file)
@@ -86,6 +86,28 @@ gmx_moltype_t::~gmx_moltype_t()
     done_atom(&atoms);
 }
 
+static int gmx_mtop_maxresnr(const gmx::ArrayRef<const gmx_moltype_t> 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<int>::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);
index a3fe344f77b20e8d0c76de4aef9d053b92e8d0e5..6ac65d0b39d104f9c8bbf9d38872eb81cc4d58b1 100644 (file)
@@ -177,10 +177,6 @@ struct gmx_mtop_t //NOLINT(clang-analyzer-optin.performance.Padding)
     std::unique_ptr<InteractionLists> 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<int> 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> 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