Move natoms_mol out of gmx_molblock_t

natoms_mol is derived information, so moved to MolecluleBlockIndices.

Change-Id: I10080dc8fd43bac1f1976551d66218fcc1d2c53c

diff --git a/src/gromacs/domdec/domdec_topology.cpp b/src/gromacs/domdec/domdec_topology.cpp
index 30a7b7d5d336f67d1f39cbcf7f7e105a62abcc0b..e8880dc9505d925a305903a3410c8977cb8d42a8 100644 (file)
--- a/src/gromacs/domdec/domdec_topology.cpp
+++ b/src/gromacs/domdec/domdec_topology.cpp
@@ -88,10 +88,12 @@
 /*! \brief The number of integer item in the local state, used for broadcasting of the state */
 #define NITEM_DD_INIT_LOCAL_STATE 5
 
-typedef struct {
-    int  *index;  /* Index for each atom into il                  */
-    int  *il;     /* ftype|type|a0|...|an|ftype|...               */
-} reverse_ilist_t;
+struct reverse_ilist_t
+{
+    int  *index;              /* Index for each atom into il          */
+    int  *il;                 /* ftype|type|a0|...|an|ftype|...       */
+    int   numAtomsInMolecule; /* The number of atoms in this molecule */
+};
 
 typedef struct {
     int  a_start;
@@ -320,13 +322,14 @@ static void print_missing_interactions_atoms(FILE *fplog, t_commrec *cr,
     int a_end = 0;
     for (const gmx_molblock_t &molb :  mtop->molblock)
     {
-        int a_start = a_end;
-        a_end       = a_start + molb.nmol*molb.natoms_mol;
+        const gmx_moltype_t &moltype  = mtop->moltype[molb.type];
+        int                  a_start  = a_end;
+        a_end                        = a_start + molb.nmol*moltype.atoms.nr;
 
         print_missing_interactions_mb(fplog, cr, rt,
-                                      *(mtop->moltype[molb.type].name),
+                                      *(moltype.name),
                                       &rt->ril_mt[molb.type],
-                                      a_start, a_end, molb.natoms_mol,
+                                      a_start, a_end, moltype.atoms.nr,
                                       molb.nmol,
                                       idef);
     }
@@ -645,6 +648,8 @@ static int make_reverse_ilist(const t_ilist *ilist,
 
     sfree(count);
 
+    ril_mt->numAtomsInMolecule = atoms->nr;
+
     return nint_mt;
 }
 
@@ -739,11 +744,13 @@ static gmx_reverse_top_t *make_reverse_top(const gmx_mtop_t *mtop, gmx_bool bFE,
     int i         = 0;
     for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
     {
-        rt->mbi[mb].a_start    = i;
-        i += mtop->molblock[mb].nmol*mtop->molblock[mb].natoms_mol;
-        rt->mbi[mb].a_end      = i;
-        rt->mbi[mb].natoms_mol = mtop->molblock[mb].natoms_mol;
-        rt->mbi[mb].type       = mtop->molblock[mb].type;
+        const gmx_molblock_t &molb           = mtop->molblock[mb];
+        int                   numAtomsPerMol = mtop->moltype[molb.type].atoms.nr;
+        rt->mbi[mb].a_start                  = i;
+        i                                   += molb.nmol*numAtomsPerMol;
+        rt->mbi[mb].a_end                    = i;
+        rt->mbi[mb].natoms_mol               = numAtomsPerMol;
+        rt->mbi[mb].type                     = molb.type;
     }
 
     rt->nthread = gmx_omp_nthreads_get(emntDomdec);
@@ -916,7 +923,8 @@ static inline void add_ifunc(int nral, t_iatom *tiatoms, t_ilist *il)
 }
 
 /*! \brief Store a position restraint in idef and iatoms, complex because the parameters are different for each entry */
-static void add_posres(int mol, int a_mol, const gmx_molblock_t *molb,
+static void add_posres(int mol, int a_mol, int numAtomsInMolecule,
+                       const gmx_molblock_t *molb,
                        t_iatom *iatoms, const t_iparams *ip_in,
                        t_idef *idef)
 {
@@ -937,11 +945,8 @@ static void add_posres(int mol, int a_mol, const gmx_molblock_t *molb,
     *ip = ip_in[iatoms[0]];
 
     /* Get the position restraint coordinates from the molblock */
-    a_molb = mol*molb->natoms_mol + a_mol;
-    if (a_molb >= static_cast<int>(molb->posres_xA.size()))
-    {
-        gmx_incons("Not enough position restraint coordinates");
-    }
+    a_molb = mol*numAtomsInMolecule + a_mol;
+    GMX_ASSERT(a_molb < static_cast<int>(molb->posres_xA.size()), "We need a sufficient number of 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];
@@ -962,7 +967,8 @@ static void add_posres(int mol, int a_mol, const gmx_molblock_t *molb,
 }
 
 /*! \brief Store a flat-bottomed position restraint in idef and iatoms, complex because the parameters are different for each entry */
-static void add_fbposres(int mol, int a_mol, const gmx_molblock_t *molb,
+static void add_fbposres(int mol, int a_mol, int numAtomsInMolecule,
+                         const gmx_molblock_t *molb,
                          t_iatom *iatoms, const t_iparams *ip_in,
                          t_idef *idef)
 {
@@ -982,12 +988,9 @@ static void add_fbposres(int mol, int a_mol, const gmx_molblock_t *molb,
     /* Copy the force constants */
     *ip = ip_in[iatoms[0]];
 
-    /* Get the position restriant coordinats from the molblock */
-    a_molb = mol*molb->natoms_mol + a_mol;
-    if (a_molb >= static_cast<int>(molb->posres_xA.size()))
-    {
-        gmx_incons("Not enough position restraint coordinates");
-    }
+    /* Get the position restraint coordinats from the molblock */
+    a_molb = mol*numAtomsInMolecule + a_mol;
+    GMX_ASSERT(a_molb < static_cast<int>(molb->posres_xA.size()), "We need a sufficient number of position restraint coordinates");
     /* Take reference positions from A position of normal posres */
     ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
     ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
@@ -1279,6 +1282,7 @@ static void combine_idef(t_idef *dest, const thread_work_t *src, int nsrc,
 static inline void
 check_assign_interactions_atom(int i, int i_gl,
                                int mol, int i_mol,
+                               int numAtomsInMolecule,
                                const int *index, const int *rtil,
                                gmx_bool bInterMolInteractions,
                                int ind_start, int ind_end,
@@ -1358,13 +1362,13 @@ check_assign_interactions_atom(int i, int i_gl,
                     tiatoms[1] = i;
                     if (ftype == F_POSRES)
                     {
-                        add_posres(mol, i_mol, molb, tiatoms, ip_in,
-                                   idef);
+                        add_posres(mol, i_mol, numAtomsInMolecule,
+                                   molb, tiatoms, ip_in, idef);
                     }
                     else if (ftype == F_FBPOSRES)
                     {
-                        add_fbposres(mol, i_mol, molb, tiatoms, ip_in,
-                                     idef);
+                        add_fbposres(mol, i_mol, numAtomsInMolecule,
+                                     molb, tiatoms, ip_in, idef);
                     }
                 }
                 else
@@ -1557,6 +1561,7 @@ static int make_bondeds_zone(gmx_domdec_t *dd,
         rtil  = rt->ril_mt[mt].il;
 
         check_assign_interactions_atom(i, i_gl, mol, i_mol,
+                                       rt->ril_mt[mt].numAtomsInMolecule,
                                        index, rtil, FALSE,
                                        index[i_mol], index[i_mol+1],
                                        dd, zones,
@@ -1579,6 +1584,7 @@ static int make_bondeds_zone(gmx_domdec_t *dd,
             rtil  = rt->ril_intermol.il;
 
             check_assign_interactions_atom(i, i_gl, mol, i_mol,
+                                           rt->ril_mt[mt].numAtomsInMolecule,
                                            index, rtil, TRUE,
                                            index[i_gl], index[i_gl + 1],
                                            dd, zones,

diff --git a/src/gromacs/fileio/tpxio.cpp b/src/gromacs/fileio/tpxio.cpp
index 61fc50b8f1148bb039559e8f40fd4fb8165f7071..89c90e38a070dd1eedcca99fbd3b05d5773348c0 100644 (file)
--- a/src/gromacs/fileio/tpxio.cpp
+++ b/src/gromacs/fileio/tpxio.cpp
@@ -2418,11 +2418,17 @@ static void do_moltype(t_fileio *fio, gmx_moltype_t *molt, gmx_bool bRead,
     do_blocka(fio, &molt->excls, bRead);
 }
 
-static void do_molblock(t_fileio *fio, gmx_molblock_t *molb, gmx_bool bRead)
+static void do_molblock(t_fileio *fio, gmx_molblock_t *molb,
+                        int numAtomsPerMolecule,
+                        gmx_bool bRead)
 {
     gmx_fio_do_int(fio, molb->type);
     gmx_fio_do_int(fio, molb->nmol);
-    gmx_fio_do_int(fio, molb->natoms_mol);
+    /* To maintain forward topology reading compatibility, we store #atoms.
+     * TODO: Change this to conditional reading of a dummy int when we
+     *       increase tpx_generation.
+     */
+    gmx_fio_do_int(fio, numAtomsPerMolecule);
     /* Position restraint coordinates */
     int numPosres_xA = molb->posres_xA.size();
     gmx_fio_do_int(fio, numPosres_xA);
@@ -2511,7 +2517,8 @@ static void do_mtop(t_fileio *fio, gmx_mtop_t *mtop, gmx_bool bRead,
     }
     for (gmx_molblock_t &molblock : mtop->molblock)
     {
-        do_molblock(fio, &molblock, bRead);
+        int numAtomsPerMolecule = (bRead ? 0 : mtop->moltype[molblock.type].atoms.nr);
+        do_molblock(fio, &molblock, numAtomsPerMolecule, bRead);
     }
     gmx_fio_do_int(fio, mtop->natoms);
 

diff --git a/src/gromacs/gmxlib/chargegroup.cpp b/src/gromacs/gmxlib/chargegroup.cpp
index 19f70e9eba9d9da2cc000afb4446de4c10eae983..fc45a1dc723b67e701bd6788d60c12d048bb0393 100644 (file)
--- a/src/gromacs/gmxlib/chargegroup.cpp
+++ b/src/gromacs/gmxlib/chargegroup.cpp
@@ -128,7 +128,7 @@ void calc_chargegroup_radii(const gmx_mtop_t *mtop, rvec *x,
                     }
                 }
             }
-            a_mol += molb.natoms_mol;
+            a_mol += molt->atoms.nr;
         }
     }
 

diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp
index d124df0596ec239f2cd7ae0078b02d449e72f2cd..5979bf5b26807244f45a572581869b369050216f 100644 (file)
--- a/src/gromacs/gmxpreprocess/grompp.cpp
+++ b/src/gromacs/gmxpreprocess/grompp.cpp
@@ -532,15 +532,14 @@ new_status(const char *topfile, const char *topppfile, const char *confin,
         {
             /* Merge consecutive blocks with the same molecule type */
             sys->molblock.back().nmol += molb.nmol;
-            sys->natoms += molb.nmol*sys->molblock.back().natoms_mol;
+            sys->natoms += molb.nmol*sys->moltype[sys->molblock.back().type].atoms.nr;
         }
         else if (molb.nmol > 0)
         {
             /* Add a new molblock to the topology */
             sys->molblock.push_back(molb);
             gmx_molblock_t &molbs  = sys->molblock.back();
-            molbs.natoms_mol       = molinfo[molbs.type].atoms.nr;
-            sys->natoms           += molbs.nmol*molbs.natoms_mol;
+            sys->natoms           += molbs.nmol*molinfo[molbs.type].atoms.nr;
         }
     }
     if (sys->molblock.empty())

diff --git a/src/gromacs/mdlib/broadcaststructs.cpp b/src/gromacs/mdlib/broadcaststructs.cpp
index bba22c6e378ae1d376304696a0a84d860e9bff43..9c70fbc1fb21e032334a3caf4268bc7d022ff397 100644 (file)
--- a/src/gromacs/mdlib/broadcaststructs.cpp
+++ b/src/gromacs/mdlib/broadcaststructs.cpp
@@ -769,7 +769,6 @@ static void bc_molblock(const t_commrec *cr, gmx_molblock_t *molb)
     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");

diff --git a/src/gromacs/mdlib/membed.cpp b/src/gromacs/mdlib/membed.cpp
index 3779c96ea85af7cef6bff78c4d57f2ec658d1946..841aa6053d824588eadbf8b1817128badfab6dfe 100644 (file)
--- a/src/gromacs/mdlib/membed.cpp
+++ b/src/gromacs/mdlib/membed.cpp
@@ -443,7 +443,7 @@ static int init_mem_at(mem_t *mem_p, gmx_mtop_t *mtop, rvec *r, matrix box, pos_
     mem_p->zmed = (zmax-zmin)/2+zmin;
 
     /*number of membrane molecules in protein box*/
-    nmolbox = count/mtop->molblock[block].natoms_mol;
+    nmolbox = count/mtop->moltype[mtop->molblock[block].type].atoms.nr;
     /*membrane area within the box defined by the min and max coordinates of the embedded molecule*/
     mem_area = (pos_ins->xmax[XX]-pos_ins->xmin[XX])*(pos_ins->xmax[YY]-pos_ins->xmin[YY]);
     /*rough estimate of area per lipid, assuming there is only one type of lipid in the membrane*/
@@ -583,7 +583,7 @@ static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc
                             {
                                 z_lip += r[k][ZZ];
                             }
-                            z_lip /= mtop->molblock[block].natoms_mol;
+                            z_lip /=  molecules.index[mol_id+1] - molecules.index[mol_id];
                             if (z_lip < mem_p->zmed)
                             {
                                 nlower++;
@@ -653,7 +653,7 @@ static int gen_rm_list(rm_t *rm_p, t_block *ins_at, t_block *rest_at, t_pbc *pbc
                 {
                     z_lip += r[k][ZZ];
                 }
-                z_lip /= mtop->molblock[block].natoms_mol;
+                z_lip /= molecules.index[mol_id+1] - molecules.index[mol_id];
                 if (nupper > nlower && z_lip < mem_p->zmed)
                 {
                     rm_p->mol[nrm]   = mol_id;
@@ -709,7 +709,7 @@ static void rm_group(gmx_groups_t *groups, gmx_mtop_t *mtop, rm_t *rm_p, t_state
         at     = molecules.index[mol_id];
         block  = rm_p->block[i];
         mtop->molblock[block].nmol--;
-        for (j = 0; j < mtop->molblock[block].natoms_mol; j++)
+        for (j = 0; j < mtop->moltype[mtop->molblock[block].type].atoms.nr; j++)
         {
             list[n] = at+j;
             n++;
@@ -832,7 +832,7 @@ static int rm_bonded(t_block *ins_at, gmx_mtop_t *mtop)
     {
         /*loop over molecule blocks*/
         type         = mtop->molblock[i].type;
-        natom        = mtop->molblock[i].natoms_mol;
+        natom        = mtop->moltype[type].atoms.nr;
         nmol         = mtop->molblock[i].nmol;
 
         for (j = 0; j < natom*nmol && bRM[type] == TRUE; j++)

diff --git a/src/gromacs/mdlib/sim_util.cpp b/src/gromacs/mdlib/sim_util.cpp
index 9d602459b0e283f8bc5d7e1a34dde256e5eaf49b..137e5028c84ff22af597968e1dac9b41ba0f968b 100644 (file)
--- a/src/gromacs/mdlib/sim_util.cpp
+++ b/src/gromacs/mdlib/sim_util.cpp
@@ -2593,17 +2593,18 @@ static void low_do_pbc_mtop(FILE *fplog, int ePBC, const matrix box,
     as = 0;
     for (const gmx_molblock_t &molb : mtop->molblock)
     {
-        if (molb.natoms_mol == 1 ||
-            (!bFirst && mtop->moltype[molb.type].cgs.nr == 1))
+        const gmx_moltype_t &moltype = mtop->moltype[molb.type];
+        if (moltype.atoms.nr == 1 ||
+            (!bFirst && moltype.cgs.nr == 1))
         {
             /* Just one atom or charge group in the molecule, no PBC required */
-            as += molb.nmol*molb.natoms_mol;
+            as += molb.nmol*moltype.atoms.nr;
         }
         else
         {
             /* Pass NULL iso fplog to avoid graph prints for each molecule type */
-            mk_graph_ilist(nullptr, mtop->moltype[molb.type].ilist,
-                           0, molb.natoms_mol, FALSE, FALSE, graph);
+            mk_graph_ilist(nullptr, moltype.ilist,
+                           0, moltype.atoms.nr, FALSE, FALSE, graph);
 
             for (mol = 0; mol < molb.nmol; mol++)
             {
@@ -2615,7 +2616,7 @@ static void low_do_pbc_mtop(FILE *fplog, int ePBC, const matrix box,
                  * since we no longer need this graph.
                  */
 
-                as += molb.natoms_mol;
+                as += moltype.atoms.nr;
             }
             done_graph(graph);
         }

diff --git a/src/gromacs/selection/indexutil.cpp b/src/gromacs/selection/indexutil.cpp
index 959a8623a4bd0d3b0e7a787659d603047b84ac62..b526ad90ebfbf5457a11a7c6c51471b3789c3759 100644 (file)
--- a/src/gromacs/selection/indexutil.cpp
+++ b/src/gromacs/selection/indexutil.cpp
@@ -942,7 +942,7 @@ gmx_ana_index_make_block(t_blocka *t, const gmx_mtop_t *top, gmx_ana_index_t *g,
                             --first_atom;
                         }
                         int first_mol_atom = top->moleculeBlockIndices[molb].globalAtomStart;
-                        first_mol_atom += molnr*top->molblock[molb].natoms_mol;
+                        first_mol_atom += molnr*top->moleculeBlockIndices[molb].numAtomsPerMolecule;
                         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)
@@ -958,10 +958,9 @@ gmx_ana_index_make_block(t_blocka *t, const gmx_mtop_t *top, gmx_ana_index_t *g,
                         {
                             ++molb;
                         }
-                        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)
+                        const int                   atomStart     = blockIndices.globalAtomStart + (id - blockIndices.moleculeIndexStart)*blockIndices.numAtomsPerMolecule;
+                        for (int j = 0; j < blockIndices.numAtomsPerMolecule; ++j)
                         {
                             t->a[t->nra++] = atomStart + j;
                         }

diff --git a/src/gromacs/selection/tests/toputils.cpp b/src/gromacs/selection/tests/toputils.cpp
index c33343e66a9c304ff17058e17aaa521de503a548..677cfb5e31d37f25263ec7dcd6c2c0ad6f16935a 100644 (file)
--- a/src/gromacs/selection/tests/toputils.cpp
+++ b/src/gromacs/selection/tests/toputils.cpp
@@ -158,7 +158,6 @@ void TopologyManager::initAtoms(int count)
     mtop_->molblock.resize(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_.get());
@@ -237,8 +236,8 @@ void TopologyManager::initUniformMolecules(int moleculeSize)
     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;
+    gmx_mtop_finalize(mtop_.get());
 }
 
 void TopologyManager::initFrameIndices(const ArrayRef<const int> &index)

diff --git a/src/gromacs/swap/swapcoords.cpp b/src/gromacs/swap/swapcoords.cpp
index 801a0574487f0db0f3b6959f4891901332109559..186e8fbf24cec01623a5470a788563ac69fd70c8 100644 (file)
--- a/src/gromacs/swap/swapcoords.cpp
+++ b/src/gromacs/swap/swapcoords.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.
@@ -989,7 +989,7 @@ static int get_group_apm_check(
      * first solvent atom: */
     int molb = 0;
     mtopGetMolblockIndex(mtop, ind[0], &molb, nullptr, nullptr);
-    int apm = mtop->molblock[molb].natoms_mol;
+    int apm = mtop->moleculeBlockIndices[molb].numAtomsPerMolecule;
 
     if (bVerbose)
     {
@@ -1001,7 +1001,7 @@ static int get_group_apm_check(
     for (int i = 1; i < nat; i++)
     {
         mtopGetMolblockIndex(mtop, ind[i], &molb, nullptr, nullptr);
-        if (apm != mtop->molblock[molb].natoms_mol)
+        if (apm != mtop->moleculeBlockIndices[molb].numAtomsPerMolecule)
         {
             gmx_fatal(FARGS, "Not all molecules of swap group %d consist of %d atoms.",
                       igroup, apm);

diff --git a/src/gromacs/tools/convert_tpr.cpp b/src/gromacs/tools/convert_tpr.cpp
index 8b5f9cbaaa880c21065a34c5459d75b731584c25..242a1a76cca1a8ce6e5f4094a287348e8e1b6f74 100644 (file)
--- a/src/gromacs/tools/convert_tpr.cpp
+++ b/src/gromacs/tools/convert_tpr.cpp
@@ -300,11 +300,10 @@ static void reduce_topology_x(int gnx, int index[],
     mtop->moltype[0].excls = top.excls;
 
     mtop->molblock.resize(1);
-    mtop->molblock[0].type       = 0;
-    mtop->molblock[0].nmol       = 1;
-    mtop->molblock[0].natoms_mol = top.atoms.nr;
+    mtop->molblock[0].type = 0;
+    mtop->molblock[0].nmol = 1;
 
-    mtop->natoms                 = top.atoms.nr;
+    mtop->natoms           = top.atoms.nr;
 }
 
 static void zeroq(int index[], gmx_mtop_t *mtop)

diff --git a/src/gromacs/topology/mtop_lookup.h b/src/gromacs/topology/mtop_lookup.h
index 1309b265d8fbecc8983ba9dc6a020bb03386dde1..d66ee55f9cb56395b386dcb180de46a7aead9e34 100644 (file)
--- a/src/gromacs/topology/mtop_lookup.h
+++ b/src/gromacs/topology/mtop_lookup.h
@@ -74,8 +74,7 @@ mtopGetMolblockIndex(const gmx_mtop_t *mtop,
 {
     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");
+    GMX_ASSERT(*moleculeBlock >= 0 && *moleculeBlock < static_cast<int>(mtop->moleculeBlockIndices.size()), "The starting molecule block index for the search should be within range and moleculeBlockIndices should not be empty");
 
     /* Search the molecue block index using bisection */
     int molBlock0 = -1;
@@ -100,14 +99,14 @@ mtopGetMolblockIndex(const gmx_mtop_t *mtop,
         *moleculeBlock = ((molBlock0 + molBlock1 + 1) >> 1);
     }
 
-    int molIndex = (globalAtomIndex - globalAtomStart) / mtop->molblock[*moleculeBlock].natoms_mol;
+    int molIndex = (globalAtomIndex - globalAtomStart) / mtop->moleculeBlockIndices[*moleculeBlock].numAtomsPerMolecule;
     if (moleculeIndex != nullptr)
     {
         *moleculeIndex = molIndex;
     }
     if (atomIndexInMolecule != nullptr)
     {
-        *atomIndexInMolecule = globalAtomIndex - globalAtomStart - molIndex*mtop->molblock[*moleculeBlock].natoms_mol;
+        *atomIndexInMolecule = globalAtomIndex - globalAtomStart - molIndex*mtop->moleculeBlockIndices[*moleculeBlock].numAtomsPerMolecule;
     }
 }
 

diff --git a/src/gromacs/topology/mtop_util.cpp b/src/gromacs/topology/mtop_util.cpp
index 4fb4fd0d1c09f6d257a235cd0dbd554ec3bdcc9a..120677ed8dfed4f27f965525b94f3cf1525651dd 100644 (file)
--- a/src/gromacs/topology/mtop_util.cpp
+++ b/src/gromacs/topology/mtop_util.cpp
@@ -90,9 +90,10 @@ static void buildMolblockIndices(gmx_mtop_t *mtop)
         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*molb.natoms_mol;
+        atomIndex                    += molb.nmol*indices.numAtomsPerMolecule;
         residueIndex                 += molb.nmol*numResPerMol;
         indices.globalAtomEnd         = atomIndex;
         indices.residueNumberStart    = residueNumberStart;
@@ -460,7 +461,7 @@ gmx_bool gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop,
 
     if (iloop->mol >= 0)
     {
-        iloop->a_offset += iloop->mtop->molblock[iloop->mblock].natoms_mol;
+        iloop->a_offset += iloop->mtop->moleculeBlockIndices[iloop->mblock].numAtomsPerMolecule;
     }
 
     iloop->mol++;
@@ -976,9 +977,10 @@ static void fillMoleculeIndices(const gmx_mtop_t  &mtop,
     index[globalMolIndex] = globalAtomIndex;
     for (const gmx_molblock_t &molb : mtop.molblock)
     {
+        int numAtomsPerMolecule = mtop.moltype[molb.type].atoms.nr;
         for (int mol = 0; mol < molb.nmol; mol++)
         {
-            globalAtomIndex       += molb.natoms_mol;
+            globalAtomIndex       += numAtomsPerMolecule;
             globalMolIndex        += 1;
             index[globalMolIndex]  = globalAtomIndex;
         }
@@ -1076,11 +1078,10 @@ void convertAtomsToMtop(t_symtab    *symtab,
     mtop->molblock.resize(1);
     mtop->molblock[0].type       = 0;
     mtop->molblock[0].nmol       = 1;
-    mtop->molblock[0].natoms_mol = mtop->moltype[mtop->molblock[0].type].atoms.nr;
 
     mtop->bIntermolecularInteractions = FALSE;
 
-    mtop->natoms                 = mtop->molblock[0].natoms_mol;
+    mtop->natoms                 = atoms->nr;
 
     mtop->haveMoleculeIndices    = false;
 

diff --git a/src/gromacs/topology/topology.cpp b/src/gromacs/topology/topology.cpp
index 591ac78ba0e747032928ddf288adbfd668b22fb1..08d6b0676e6a8150b674ce06eea8ba7e07cac0d8 100644 (file)
--- a/src/gromacs/topology/topology.cpp
+++ b/src/gromacs/topology/topology.cpp
@@ -361,7 +361,6 @@ static void pr_molblock(FILE *fp, int indent, const char *title,
     fprintf(fp, "%-20s = %d \"%s\"\n",
             "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->posres_xA.size());
     if (!molb->posres_xA.empty())
     {

diff --git a/src/gromacs/topology/topology.h b/src/gromacs/topology/topology.h
index 7d9affdf2c4068087b40be58ccfe30c787426c0b..fd35fb4bad0176c18c6610bb26c0636cf9735c4b 100644 (file)
--- a/src/gromacs/topology/topology.h
+++ b/src/gromacs/topology/topology.h
@@ -85,19 +85,17 @@ struct gmx_molblock_t
     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 */
-
-    /* Convenience information, derived from other gmx_mtop_t contents */
-    int     natoms_mol = 0;           /**< 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 */
-    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 */
+    int     numAtomsPerMolecule; /**< Number of atoms in a molecule in the block */
+    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 */
+    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 */
 };
 
 typedef struct gmx_groups_t