Refactor t_molinfo to MoleculeInformation class
authorPaul Bauer <paul.bauer.q@gmail.com>
Wed, 13 Feb 2019 17:46:15 +0000 (18:46 +0100)
committerMagnus Lundborg <magnus.lundborg@scilifelab.se>
Thu, 21 Feb 2019 09:55:40 +0000 (10:55 +0100)
Refs #2833

Change-Id: I66a6afc8ade636fb3002492e6c4e1d30ae51019a

12 files changed:
src/gromacs/gmxpreprocess/convparm.cpp
src/gromacs/gmxpreprocess/convparm.h
src/gromacs/gmxpreprocess/grompp.cpp
src/gromacs/gmxpreprocess/grompp_impl.h
src/gromacs/gmxpreprocess/tomorse.cpp
src/gromacs/gmxpreprocess/tomorse.h
src/gromacs/gmxpreprocess/topio.cpp
src/gromacs/gmxpreprocess/topio.h
src/gromacs/gmxpreprocess/toppush.cpp
src/gromacs/gmxpreprocess/toppush.h
src/gromacs/gmxpreprocess/toputil.cpp
src/gromacs/gmxpreprocess/toputil.h

index df6889d52a3b6c61e48df38b907a1066551cecc3..1ea48876c9bf4f220b20a5bcedb0de92d7c1a070 100644 (file)
@@ -497,7 +497,7 @@ static void append_interaction(InteractionList *ilist,
     }
 }
 
-static void enter_function(t_params *p, t_functype ftype, int comb, real reppow,
+static void enter_function(const t_params *p, t_functype ftype, int comb, real reppow,
                            gmx_ffparams_t *ffparams, InteractionList *il,
                            bool bNB, bool bAppend)
 {
@@ -520,7 +520,8 @@ static void enter_function(t_params *p, t_functype ftype, int comb, real reppow,
 }
 
 void convert_params(int atnr, t_params nbtypes[],
-                    t_molinfo *mi, t_molinfo *intermolecular_interactions,
+                    gmx::ArrayRef<const MoleculeInformation> mi,
+                    const MoleculeInformation *intermolecular_interactions,
                     int comb, double reppow, real fudgeQQ,
                     gmx_mtop_t *mtop)
 {
@@ -528,7 +529,6 @@ void convert_params(int atnr, t_params nbtypes[],
     unsigned long   flags;
     gmx_ffparams_t *ffp;
     gmx_moltype_t  *molt;
-    t_params       *plist;
 
     ffp           = &mtop->ffparams;
     ffp->atnr     = atnr;
@@ -548,7 +548,7 @@ void convert_params(int atnr, t_params nbtypes[],
         {
             molt->ilist[i].iatoms.clear();
 
-            plist = mi[mt].plist;
+            const t_params *plist = mi[mt].plist;
 
             flags = interaction_function[i].flags;
             if ((i != F_LJ) && (i != F_BHAM) && ((flags & IF_BOND) ||
@@ -572,7 +572,7 @@ void convert_params(int atnr, t_params nbtypes[],
         {
             (*mtop->intermolecular_ilist)[i].iatoms.clear();
 
-            plist = intermolecular_interactions->plist;
+            const t_params *plist = intermolecular_interactions->plist;
 
             if (plist[i].nr > 0)
             {
index 496d866d3d4eaca725a245a2be053f69b346e91e..2b61d6fcfd152844a0196df5679f36c0b4992c53 100644 (file)
 #ifndef GMX_GMXPREPROCESS_CONVPARM_H
 #define GMX_GMXPREPROCESS_CONVPARM_H
 
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/real.h"
 
 struct gmx_mtop_t;
-struct t_molinfo;
+struct MoleculeInformation;
 struct t_params;
 
 void convert_params(int atnr, t_params nbtypes[],
-                    t_molinfo *mi,
-                    t_molinfo *intermolecular_interactions,
+                    gmx::ArrayRef<const MoleculeInformation> mi,
+                    const MoleculeInformation *intermolecular_interactions,
                     int comb, double reppow, real fudgeQQ,
                     gmx_mtop_t *mtop);
 
index 08ce54c750b9b5404ba3388005bc0dbcf3901370..c7d7ffa6eb9df677e8a8db13c116d07e12192c51 100644 (file)
 #include "gromacs/utility/smalloc.h"
 #include "gromacs/utility/snprintf.h"
 
-static int rm_interactions(int ifunc, int nrmols, t_molinfo mols[])
+void MoleculeInformation::initMolInfo()
 {
-    int  i, n;
+    init_plist(plist);
+    init_block(&cgs);
+    init_block(&mols);
+    init_blocka(&excls);
+    init_t_atoms(&atoms, 0, FALSE);
+}
+
+void MoleculeInformation::partialCleanUp()
+{
+    done_block(&mols);
+    done_plist(plist);
+}
 
-    n = 0;
+void MoleculeInformation::fullCleanUp()
+{
+    done_atom (&atoms);
+    done_block(&cgs);
+    done_block(&mols);
+    done_plist(plist);
+}
+
+static int rm_interactions(int ifunc, gmx::ArrayRef<MoleculeInformation> mols)
+{
+    int n = 0;
     /* For all the molecule types */
-    for (i = 0; i < nrmols; i++)
+    for (auto &mol : mols)
     {
-        n += mols[i].plist[ifunc].nr;
-        mols[i].plist[ifunc].nr = 0;
+        n                  += mol.plist[ifunc].nr;
+        mol.plist[ifunc].nr = 0;
     }
     return n;
 }
@@ -406,7 +427,7 @@ static void check_shells_inputrec(gmx_mtop_t *mtop,
 
 /* TODO Decide whether this function can be consolidated with
  * gmx_mtop_ftype_count */
-static int nint_ftype(gmx_mtop_t *mtop, t_molinfo *mi, int ftype)
+static int nint_ftype(gmx_mtop_t *mtop, gmx::ArrayRef<const MoleculeInformation> mi, int ftype)
 {
     int nint = 0;
     for (const gmx_molblock_t &molb : mtop->molblock)
@@ -421,73 +442,63 @@ static int nint_ftype(gmx_mtop_t *mtop, t_molinfo *mi, int ftype)
  * in the order of use in the molblocks,
  * unused molecule types are deleted.
  */
-static void renumber_moltypes(gmx_mtop_t *sys,
-                              int *nmolinfo, t_molinfo **molinfo)
+static void renumber_moltypes(gmx_mtop_t                       *sys,
+                              std::vector<MoleculeInformation> *molinfo)
 {
-    int       *order, norder;
-    t_molinfo *minew;
 
-    snew(order, *nmolinfo);
-    norder = 0;
+    std::vector<int> order;
     for (gmx_molblock_t &molblock : sys->molblock)
     {
-        int i;
-        for (i = 0; i < norder; i++)
+        const auto found = std::find_if(order.begin(), order.end(),
+                                        [&molblock](const auto &entry)
+                                        { return molblock.type == entry; });
+        if (found == order.end())
         {
-            if (order[i] == molblock.type)
-            {
-                break;
-            }
+            /* This type did not occur yet, add it */
+            order.push_back(molblock.type);
+            molblock.type = order.size() - 1;
         }
-        if (i == norder)
+        else
         {
-            /* This type did not occur yet, add it */
-            order[norder] = molblock.type;
-            /* Renumber the moltype in the topology */
-            norder++;
+            molblock.type = std::distance(order.begin(), found);
         }
-        molblock.type = i;
     }
 
     /* We still need to reorder the molinfo structs */
-    snew(minew, norder);
-    for (int mi = 0; mi < *nmolinfo; mi++)
+    std::vector<MoleculeInformation> minew(order.size());
+    int index = 0;
+    for (auto &mi : *molinfo)
     {
-        int i;
-        for (i = 0; i < norder; i++)
+        const auto found = std::find(order.begin(), order.end(), index);
+        if (found != order.end())
         {
-            if (order[i] == mi)
-            {
-                break;
-            }
-        }
-        if (i == norder)
-        {
-            done_mi(&(*molinfo)[mi]);
+            int pos = std::distance(order.begin(), found);
+            minew[pos] = mi;
         }
         else
         {
-            minew[i] = (*molinfo)[mi];
+            // Need to manually clean up memory ....
+            mi.fullCleanUp();
         }
+        index++;
     }
-    sfree(order);
-    sfree(*molinfo);
 
-    *nmolinfo = norder;
     *molinfo  = minew;
 }
 
-static void molinfo2mtop(int nmi, t_molinfo *mi, gmx_mtop_t *mtop)
+static void molinfo2mtop(gmx::ArrayRef<const MoleculeInformation> mi, gmx_mtop_t *mtop)
 {
-    mtop->moltype.resize(nmi);
-    for (int m = 0; m < nmi; m++)
+    mtop->moltype.resize(mi.size());
+    int pos = 0;
+    for (const auto &mol : mi)
     {
-        gmx_moltype_t &molt = mtop->moltype[m];
-        molt.name           = mi[m].name;
-        molt.atoms          = mi[m].atoms;
+        gmx_moltype_t &molt = mtop->moltype[pos];
+        molt.name           = mol.name;
+        molt.atoms          = mol.atoms;
         /* ilists are copied later */
-        molt.cgs            = mi[m].cgs;
-        molt.excls          = mi[m].excls;
+        molt.cgs            = mol.cgs;
+        molt.excls          = mol.excls;
+        pos++;
     }
 }
 
@@ -496,23 +507,23 @@ new_status(const char *topfile, const char *topppfile, const char *confin,
            t_gromppopts *opts, t_inputrec *ir, gmx_bool bZero,
            bool bGenVel, bool bVerbose, t_state *state,
            gpp_atomtype *atype, gmx_mtop_t *sys,
-           int *nmi, t_molinfo **mi, t_molinfo **intermolecular_interactions,
+           std::vector<MoleculeInformation> *mi,
+           std::unique_ptr<MoleculeInformation> *intermolecular_interactions,
            t_params plist[],
            int *comb, double *reppow, real *fudgeQQ,
            gmx_bool bMorse,
            warninp *wi)
 {
-    t_molinfo                  *molinfo = nullptr;
-    std::vector<gmx_molblock_t> molblock;
-    int                         i, nrmols, nmismatch;
-    bool                        ffParametrizedWithHBondConstraints;
-    char                        buf[STRLEN];
-    char                        warn_buf[STRLEN];
+    std::vector<gmx_molblock_t>      molblock;
+    int                              i, nrmols, nmismatch;
+    bool                             ffParametrizedWithHBondConstraints;
+    char                             buf[STRLEN];
+    char                             warn_buf[STRLEN];
 
     /* TOPOLOGY processing */
     sys->name = do_top(bVerbose, topfile, topppfile, opts, bZero, &(sys->symtab),
                        plist, comb, reppow, fudgeQQ,
-                       atype, &nrmols, &molinfo, intermolecular_interactions,
+                       atype, &nrmols, mi, intermolecular_interactions,
                        ir,
                        &molblock,
                        &ffParametrizedWithHBondConstraints,
@@ -534,23 +545,23 @@ new_status(const char *topfile, const char *topppfile, const char *confin,
             /* Add a new molblock to the topology */
             sys->molblock.push_back(molb);
         }
-        sys->natoms += molb.nmol*molinfo[sys->molblock.back().type].atoms.nr;
+        sys->natoms += molb.nmol*(*mi)[sys->molblock.back().type].atoms.nr;
     }
     if (sys->molblock.empty())
     {
         gmx_fatal(FARGS, "No molecules were defined in the system");
     }
 
-    renumber_moltypes(sys, &nrmols, &molinfo);
+    renumber_moltypes(sys, mi);
 
     if (bMorse)
     {
-        convert_harmonics(nrmols, molinfo, atype);
+        convert_harmonics(*mi, atype);
     }
 
     if (ir->eDisre == edrNone)
     {
-        i = rm_interactions(F_DISRES, nrmols, molinfo);
+        i = rm_interactions(F_DISRES, *mi);
         if (i > 0)
         {
             set_warning_line(wi, "unknown", -1);
@@ -560,7 +571,7 @@ new_status(const char *topfile, const char *topppfile, const char *confin,
     }
     if (!opts->bOrire)
     {
-        i = rm_interactions(F_ORIRES, nrmols, molinfo);
+        i = rm_interactions(F_ORIRES, *mi);
         if (i > 0)
         {
             set_warning_line(wi, "unknown", -1);
@@ -570,7 +581,7 @@ new_status(const char *topfile, const char *topppfile, const char *confin,
     }
 
     /* Copy structures from msys to sys */
-    molinfo2mtop(nrmols, molinfo, sys);
+    molinfo2mtop(*mi, sys);
 
     gmx_mtop_finalize(sys);
 
@@ -654,9 +665,9 @@ new_status(const char *topfile, const char *topppfile, const char *confin,
         fprintf(stderr, "double-checking input for internal consistency...\n");
     }
     {
-        bool bHasNormalConstraints = 0 < (nint_ftype(sys, molinfo, F_CONSTR) +
-                                          nint_ftype(sys, molinfo, F_CONSTRNC));
-        bool bHasAnyConstraints = bHasNormalConstraints || 0 < nint_ftype(sys, molinfo, F_SETTLE);
+        bool bHasNormalConstraints = 0 < (nint_ftype(sys, *mi, F_CONSTR) +
+                                          nint_ftype(sys, *mi, F_CONSTRNC));
+        bool bHasAnyConstraints = bHasNormalConstraints || 0 < nint_ftype(sys, *mi, F_SETTLE);
         double_check(ir, state->box,
                      bHasNormalConstraints,
                      bHasAnyConstraints,
@@ -686,9 +697,6 @@ new_status(const char *topfile, const char *topppfile, const char *confin,
         stop_cm(stdout, state->natoms, mass, state->x.rvec_array(), state->v.rvec_array());
         sfree(mass);
     }
-
-    *nmi = nrmols;
-    *mi  = molinfo;
 }
 
 static void copy_state(const char *slog, t_trxframe *fr,
@@ -811,7 +819,9 @@ static void cont_status(const char *slog, const char *ener,
     }
 }
 
-static void read_posres(gmx_mtop_t *mtop, t_molinfo *molinfo, gmx_bool bTopB,
+static void read_posres(gmx_mtop_t *mtop,
+                        gmx::ArrayRef<const MoleculeInformation> molinfo,
+                        gmx_bool bTopB,
                         const char *fn,
                         int rc_scaling, int ePBC,
                         rvec com,
@@ -826,7 +836,6 @@ static void read_posres(gmx_mtop_t *mtop, t_molinfo *molinfo, gmx_bool bTopB,
     int                 natoms, npbcdim = 0;
     char                warn_buf[STRLEN];
     int                 a, i, ai, j, k, nat_molb;
-    t_params           *pr, *prfb;
     t_atom             *atom;
 
     snew(top, 1);
@@ -861,8 +870,8 @@ static void read_posres(gmx_mtop_t *mtop, t_molinfo *molinfo, gmx_bool bTopB,
     for (gmx_molblock_t &molb : mtop->molblock)
     {
         nat_molb = molb.nmol*mtop->moltype[molb.type].atoms.nr;
-        pr       = &(molinfo[molb.type].plist[F_POSRES]);
-        prfb     = &(molinfo[molb.type].plist[F_FBPOSRES]);
+        const t_params *pr   = &(molinfo[molb.type].plist[F_POSRES]);
+        const t_params *prfb = &(molinfo[molb.type].plist[F_FBPOSRES]);
         if (pr->nr > 0 || prfb->nr > 0)
         {
             atom = mtop->moltype[molb.type].atoms.atom;
@@ -988,13 +997,14 @@ static void read_posres(gmx_mtop_t *mtop, t_molinfo *molinfo, gmx_bool bTopB,
     sfree(hadAtom);
 }
 
-static void gen_posres(gmx_mtop_t *mtop, t_molinfo *mi,
+static void gen_posres(gmx_mtop_t *mtop,
+                       gmx::ArrayRef<const MoleculeInformation> mi,
                        const char *fnA, const char *fnB,
                        int rc_scaling, int ePBC,
                        rvec com, rvec comB,
                        warninp *wi)
 {
-    read_posres  (mtop, mi, FALSE, fnA, rc_scaling, ePBC, com, wi);
+    read_posres(mtop, mi, FALSE, fnA, rc_scaling, ePBC, com, wi);
     /* It is safer to simply read the b-state posres rather than trying
      * to be smart and copy the positions.
      */
@@ -1022,7 +1032,7 @@ static void set_wall_atomtype(gpp_atomtype *at, t_gromppopts *opts,
     }
 }
 
-static int nrdf_internal(t_atoms *atoms)
+static int nrdf_internal(const t_atoms *atoms)
 {
     int i, nmass, nrdf;
 
@@ -1192,17 +1202,18 @@ static void init_cmap_grid(gmx_cmap_t *cmap_grid, int ngrid, int grid_spacing)
 }
 
 
-static int count_constraints(const gmx_mtop_t *mtop, t_molinfo *mi, warninp *wi)
+static int count_constraints(const gmx_mtop_t                        *mtop,
+                             gmx::ArrayRef<const MoleculeInformation> mi,
+                             warninp                                 *wi)
 {
     int             count, count_mol, i;
-    t_params       *plist;
     char            buf[STRLEN];
 
     count = 0;
     for (const gmx_molblock_t &molb : mtop->molblock)
     {
         count_mol = 0;
-        plist     = mi[molb.type].plist;
+        const t_params *plist = mi[molb.type].plist;
 
         for (i = 0; i < F_NRE; i++)
         {
@@ -1580,7 +1591,7 @@ static void set_verlet_buffer(const gmx_mtop_t *mtop,
 
 int gmx_grompp(int argc, char *argv[])
 {
-    const char            *desc[] = {
+    const char                          *desc[] = {
         "[THISMODULE] (the gromacs preprocessor)",
         "reads a molecular topology file, checks the validity of the",
         "file, expands the topology from a molecular description to an atomic",
@@ -1678,24 +1689,24 @@ int gmx_grompp(int argc, char *argv[])
         "interpret the output messages before attempting to bypass them with",
         "this option."
     };
-    t_gromppopts          *opts;
-    int                    nmi;
-    t_molinfo             *mi, *intermolecular_interactions;
-    gpp_atomtype          *atype;
-    int                    nvsite, comb;
-    t_params              *plist;
-    real                   fudgeQQ;
-    double                 reppow;
-    const char            *mdparin;
-    int                    ntype;
-    bool                   bNeedVel, bGenVel;
-    gmx_bool               have_atomnumber;
-    gmx_output_env_t      *oenv;
-    gmx_bool               bVerbose = FALSE;
-    warninp               *wi;
-    char                   warn_buf[STRLEN];
-
-    t_filenm               fnm[] = {
+    t_gromppopts                        *opts;
+    std::vector<MoleculeInformation>     mi;
+    std::unique_ptr<MoleculeInformation> intermolecular_interactions;
+    gpp_atomtype                        *atype;
+    int                                  nvsite, comb;
+    t_params                            *plist;
+    real                                 fudgeQQ;
+    double                               reppow;
+    const char                          *mdparin;
+    int                                  ntype;
+    bool                                 bNeedVel, bGenVel;
+    gmx_bool                             have_atomnumber;
+    gmx_output_env_t                    *oenv;
+    gmx_bool                             bVerbose = FALSE;
+    warninp                             *wi;
+    char                                 warn_buf[STRLEN];
+
+    t_filenm                             fnm[] = {
         { efMDP, nullptr,  nullptr,        ffREAD  },
         { efMDP, "-po", "mdout",     ffWRITE },
         { efSTX, "-c",  nullptr,        ffREAD  },
@@ -1806,7 +1817,7 @@ int gmx_grompp(int argc, char *argv[])
     t_state state;
     new_status(fn, opt2fn_null("-pp", NFILE, fnm), opt2fn("-c", NFILE, fnm),
                opts, ir, bZero, bGenVel, bVerbose, &state,
-               atype, &sys, &nmi, &mi, &intermolecular_interactions,
+               atype, &sys, &mi, &intermolecular_interactions,
                plist, &comb, &reppow, &fudgeQQ,
                opts->bMorse,
                wi);
@@ -1842,7 +1853,7 @@ int gmx_grompp(int argc, char *argv[])
         gmx_mtop_remove_chargegroups(&sys);
     }
 
-    if (count_constraints(&sys, mi, wi) && (ir->eConstrAlg == econtSHAKE))
+    if ((count_constraints(&sys, mi, wi) != 0) && (ir->eConstrAlg == econtSHAKE))
     {
         if (ir->eI == eiCG || ir->eI == eiLBFGS)
         {
@@ -1976,7 +1987,7 @@ int gmx_grompp(int argc, char *argv[])
     }
 
     ntype = get_atomtype_ntypes(atype);
-    convert_params(ntype, plist, mi, intermolecular_interactions,
+    convert_params(ntype, plist, mi, intermolecular_interactions.get(),
                    comb, reppow, fudgeQQ, &sys);
 
     if (debug)
@@ -2314,14 +2325,12 @@ int gmx_grompp(int argc, char *argv[])
     sfree(opts);
     done_plist(plist);
     sfree(plist);
-    for (int i = 0; i < nmi; i++)
+    for (auto &mol : mi)
     {
         // Some of the contents of molinfo have been stolen, so
-        // done_mi can't be called.
-        done_block(&mi[i].mols);
-        done_plist(mi[i].plist);
+        // fullCleanUp can't be called.
+        mol.partialCleanUp();
     }
-    sfree(mi);
     done_atomtype(atype);
     done_inputrec_strings();
     output_env_done(oenv);
index 62470454d94808767d1437c3f93e4ff36ace8633..8c6489a15ff1221e3137e2a4d34497fcce57625c 100644 (file)
@@ -94,17 +94,56 @@ struct t_excls
     int           *e;       /* The excluded atoms                   */
 };
 
-struct t_molinfo
+
+/*! \libinternal \brief
+ * Holds the molecule information during preprocessing.
+ */
+struct MoleculeInformation
 {
-    char            **name;
-    int               nrexcl;       /* Number of exclusions per atom   */
-    bool              excl_set;     /* Have exclusions been generated? */
-    bool              bProcessed;   /* Has the mol been processed           */
-    t_atoms           atoms;        /* Atoms                                */
-    t_block           cgs;          /* Charge groups                        */
-    t_block           mols;         /* Molecules                            */
-    t_blocka          excls;        /* Exclusions                           */
-    t_params          plist[F_NRE]; /* Parameters in old style              */
+    //! Name of the molecule.
+    char            **name = nullptr;
+    //!Number of exclusions per atom.
+    int               nrexcl = 0;
+    //! Have exclusions been generated?.
+    bool              excl_set = false;
+    //! Has the mol been processed.
+    bool              bProcessed = false;
+    //! Atoms in the moelcule.
+    t_atoms           atoms;
+    //! Charge groups in the molecule
+    t_block           cgs;
+    //! Molecules separated in datastructure.
+    t_block           mols;
+    //! Exclusions in the molecule.
+    t_blocka          excls;
+    //! Parameters in old style.
+    t_params          plist[F_NRE];
+
+    /*! \brief
+     * Initializer.
+     *
+     * This should be removed as soon as the underlying datastructures
+     * have been cleaned up to use proper initialization and can be copy
+     * constructed.
+     */
+    void initMolInfo();
+
+    /*! \brief
+     * Partial clean up function.
+     *
+     * Should be removed once this datastructure actually owns all its own memory and
+     * elements of it are not stolen by other structures and properly copy constructed
+     * or moved.
+     * Cleans up the mols and plist datastructures but not cgs and excls.
+     */
+    void partialCleanUp();
+
+    /*! \brief
+     * Full clean up function.
+     *
+     * Should be removed once the destructor can always do this.
+     */
+    void fullCleanUp();
 };
 
 struct t_mols
index 90b296274648f90130d60b32b0691b95c3f5ef0d..9bd2c0926217d01f3354a743862f1d5f1a6849ab 100644 (file)
@@ -187,14 +187,11 @@ static real search_e_diss(int n2m, t_2morse t2m[], char *ai, char *aj)
     }
 }
 
-void convert_harmonics(int nrmols, t_molinfo mols[], gpp_atomtype *atype)
+void convert_harmonics(gmx::ArrayRef<MoleculeInformation> mols, gpp_atomtype *atype)
 {
     int       n2m;
     t_2morse *t2m;
 
-    int       i, j, k, last, ni, nj;
-    int       nrharm, nrmorse, bb;
-    real      edis, kb, b0, beta;
     bool     *bRemoveHarm;
 
     /* First get the data */
@@ -206,61 +203,63 @@ void convert_harmonics(int nrmols, t_molinfo mols[], gpp_atomtype *atype)
     }
 
     /* For all the molecule types */
-    for (i = 0; (i < nrmols); i++)
+    int i = 0;
+    for (auto &mol : mols)
     {
         /* Check how many morse and harmonic BONDSs there are, increase size of
          * morse with the number of harmonics
          */
-        nrmorse = mols[i].plist[F_MORSE].nr;
+        int nrmorse = mol.plist[F_MORSE].nr;
 
-        for (bb = 0; (bb < F_NRE); bb++)
+        for (int bb = 0; (bb < F_NRE); bb++)
         {
             if ((interaction_function[bb].flags & IF_BTYPE) && (bb != F_MORSE))
             {
-                nrharm  = mols[i].plist[bb].nr;
-                pr_alloc(nrharm, &(mols[i].plist[F_MORSE]));
+                int nrharm  = mol.plist[bb].nr;
+                pr_alloc(nrharm, &(mol.plist[F_MORSE]));
                 snew(bRemoveHarm, nrharm);
 
                 /* Now loop over the harmonics, trying to convert them */
-                for (j = 0; (j < nrharm); j++)
+                for (int j = 0; (j < nrharm); j++)
                 {
-                    ni   = mols[i].plist[bb].param[j].ai();
-                    nj   = mols[i].plist[bb].param[j].aj();
-                    edis =
+                    int  ni   = mol.plist[bb].param[j].ai();
+                    int  nj   = mol.plist[bb].param[j].aj();
+                    real edis =
                         search_e_diss(n2m, t2m,
-                                      get_atomtype_name(mols[i].atoms.atom[ni].type, atype),
-                                      get_atomtype_name(mols[i].atoms.atom[nj].type, atype));
+                                      get_atomtype_name(mol.atoms.atom[ni].type, atype),
+                                      get_atomtype_name(mol.atoms.atom[nj].type, atype));
                     if (edis != 0)
                     {
-                        bRemoveHarm[j] = TRUE;
-                        b0             = mols[i].plist[bb].param[j].c[0];
-                        kb             = mols[i].plist[bb].param[j].c[1];
-                        beta           = std::sqrt(kb/(2*edis));
-                        mols[i].plist[F_MORSE].param[nrmorse].a[0] = ni;
-                        mols[i].plist[F_MORSE].param[nrmorse].a[1] = nj;
-                        mols[i].plist[F_MORSE].param[nrmorse].c[0] = b0;
-                        mols[i].plist[F_MORSE].param[nrmorse].c[1] = edis;
-                        mols[i].plist[F_MORSE].param[nrmorse].c[2] = beta;
+                        bRemoveHarm[j] = true;
+                        real b0             = mol.plist[bb].param[j].c[0];
+                        real kb             = mol.plist[bb].param[j].c[1];
+                        real beta           = std::sqrt(kb/(2*edis));
+                        mol.plist[F_MORSE].param[nrmorse].a[0] = ni;
+                        mol.plist[F_MORSE].param[nrmorse].a[1] = nj;
+                        mol.plist[F_MORSE].param[nrmorse].c[0] = b0;
+                        mol.plist[F_MORSE].param[nrmorse].c[1] = edis;
+                        mol.plist[F_MORSE].param[nrmorse].c[2] = beta;
                         nrmorse++;
                     }
                 }
-                mols[i].plist[F_MORSE].nr = nrmorse;
+                mol.plist[F_MORSE].nr = nrmorse;
 
                 /* Now remove the harmonics */
-                for (j = last = 0; (j < nrharm); j++)
+                int last = 0;
+                for (int j = 0; (j < nrharm); j++)
                 {
                     if (!bRemoveHarm[j])
                     {
                         /* Copy it to the last position */
-                        for (k = 0; (k < MAXATOMLIST); k++)
+                        for (int k = 0; (k < MAXATOMLIST); k++)
                         {
-                            mols[i].plist[bb].param[last].a[k] =
-                                mols[i].plist[bb].param[j].a[k];
+                            mol.plist[bb].param[last].a[k] =
+                                mol.plist[bb].param[j].a[k];
                         }
-                        for (k = 0; (k < MAXFORCEPARAM); k++)
+                        for (int k = 0; (k < MAXFORCEPARAM); k++)
                         {
-                            mols[i].plist[bb].param[last].c[k] =
-                                mols[i].plist[bb].param[j].c[k];
+                            mol.plist[bb].param[last].c[k] =
+                                mol.plist[bb].param[j].c[k];
                         }
                         last++;
                     }
@@ -268,9 +267,10 @@ void convert_harmonics(int nrmols, t_molinfo mols[], gpp_atomtype *atype)
                 sfree(bRemoveHarm);
                 fprintf(stderr, "Converted %d out of %d %s to morse bonds for mol %d\n",
                         nrharm-last, nrharm, interaction_function[bb].name, i);
-                mols[i].plist[bb].nr = last;
+                mol.plist[bb].nr = last;
             }
         }
+        i++;
     }
     sfree(t2m);
 }
index 778e477cbcf27b6dcc08b236ef9a7a50f26d6c7a..9d2f510493261e4d7c5d5d2b5334350082a047af 100644 (file)
 #ifndef GMX_GMXPREPROCESS_TOMORSE_H
 #define GMX_GMXPREPROCESS_TOMORSE_H
 
+#include "gromacs/utility/arrayref.h"
+
 struct gpp_atomtype;
-struct t_molinfo;
+struct MoleculeInformation;
 
-void convert_harmonics(int nrmols, t_molinfo mols[], gpp_atomtype *atype);
+void convert_harmonics(gmx::ArrayRef<MoleculeInformation> mols, gpp_atomtype *atype);
 
 #endif
index 97a5ff7f9939ccf0128d093e31171e2a255963b3..5878613420cd4490ef30658c8ea8ab736633e842 100644 (file)
@@ -47,6 +47,7 @@
 #include <cstring>
 
 #include <algorithm>
+#include <memory>
 
 #include <unordered_set>
 #include <sys/types.h>
@@ -349,9 +350,9 @@ static char ** cpp_opts(const char *define, const char *include,
 }
 
 
-static void make_atoms_sys(const std::vector<gmx_molblock_t> &molblock,
-                           const t_molinfo                   *molinfo,
-                           t_atoms                           *atoms)
+static void make_atoms_sys(gmx::ArrayRef<const gmx_molblock_t>      molblock,
+                           gmx::ArrayRef<const MoleculeInformation> molinfo,
+                           t_atoms                                 *atoms)
 {
     atoms->nr   = 0;
     atoms->atom = nullptr;
@@ -378,8 +379,8 @@ static char **read_topol(const char *infile, const char *outfile,
                          t_symtab    *symtab,
                          gpp_atomtype *atype,
                          int         *nrmols,
-                         t_molinfo   **molinfo,
-                         t_molinfo   **intermolecular_interactions,
+                         std::vector<MoleculeInformation> *molinfo,
+                         std::unique_ptr<MoleculeInformation> *intermolecular_interactions,
                          t_params    plist[],
                          int         *combination_rule,
                          double      *reppow,
@@ -392,32 +393,32 @@ static char **read_topol(const char *infile, const char *outfile,
                          bool        usingFullRangeElectrostatics,
                          warninp    *wi)
 {
-    FILE                 *out;
-    int                   i, sl, nb_funct;
-    char                 *pline = nullptr, **title = nullptr;
-    char                  line[STRLEN], errbuf[256], comb_str[256], nb_str[256];
-    char                  genpairs[32];
-    char                 *dirstr, *dummy2;
-    int                   nrcopies, nmol, nscan, ncombs, ncopy;
-    double                fLJ, fQQ, fPOW;
-    t_molinfo            *mi0   = nullptr;
-    DirStack             *DS;
-    Directive             d, newd;
-    t_nbparam           **nbparam, **pair;
-    gmx::ExclusionBlocks *exclusionBlocks;
-    real                  fudgeLJ = -1;    /* Multiplication factor to generate 1-4 from LJ */
-    bool                  bReadDefaults, bReadMolType, bGenPairs, bWarn_copy_A_B;
-    double                qt = 0, qBt = 0; /* total charge */
-    gpp_bond_atomtype    *batype;
-    int                   lastcg = -1;
-    int                   dcatt  = -1, nmol_couple;
+    FILE                           *out;
+    int                             sl, nb_funct;
+    char                           *pline = nullptr, **title = nullptr;
+    char                            line[STRLEN], errbuf[256], comb_str[256], nb_str[256];
+    char                            genpairs[32];
+    char                           *dirstr, *dummy2;
+    int                             nrcopies, nmol, nscan, ncombs, ncopy;
+    double                          fLJ, fQQ, fPOW;
+    MoleculeInformation            *mi0   = nullptr;
+    DirStack                       *DS;
+    Directive                       d, newd;
+    t_nbparam                     **nbparam, **pair;
+    gmx::ExclusionBlocks           *exclusionBlocks;
+    real                            fudgeLJ = -1;    /* Multiplication factor to generate 1-4 from LJ */
+    bool                            bReadDefaults, bReadMolType, bGenPairs, bWarn_copy_A_B;
+    double                          qt = 0, qBt = 0; /* total charge */
+    gpp_bond_atomtype              *batype;
+    int                             lastcg = -1;
+    int                             dcatt  = -1, nmol_couple;
     /* File handling variables */
-    int                   status;
-    bool                  done;
-    gmx_cpp_t             handle;
-    char                 *tmp_line = nullptr;
-    char                  warn_buf[STRLEN];
-    const char           *floating_point_arithmetic_tip =
+    int                             status;
+    bool                            done;
+    gmx_cpp_t                       handle;
+    char                           *tmp_line = nullptr;
+    char                            warn_buf[STRLEN];
+    const char                     *floating_point_arithmetic_tip =
         "Total charge should normally be an integer. See\n"
         "http://www.gromacs.org/Documentation/Floating_Point_Arithmetic\n"
         "for discussion on how close it should be to an integer.\n";
@@ -452,8 +453,6 @@ static char **read_topol(const char *infile, const char *outfile,
 
     *reppow  = 12.0;      /* Default value for repulsion power     */
 
-    *intermolecular_interactions = nullptr;
-
     /* Init the number of CMAP torsion angles  and grid spacing */
     plist[F_CMAP].grid_spacing = 0;
     plist[F_CMAP].nc           = 0;
@@ -578,9 +577,9 @@ static char **read_topol(const char *infile, const char *outfile,
                                  * to process the intermolecular interactions
                                  * by making a "molecule" of the size of the system.
                                  */
-                                snew(*intermolecular_interactions, 1);
-                                init_molinfo(*intermolecular_interactions);
-                                mi0 = *intermolecular_interactions;
+                                *intermolecular_interactions = std::make_unique<MoleculeInformation>( );
+                                mi0 = intermolecular_interactions->get();
+                                mi0->initMolInfo();
                                 make_atoms_sys(*molblock, *molinfo,
                                                &mi0->atoms);
                             }
@@ -723,10 +722,11 @@ static char **read_topol(const char *infile, const char *outfile,
                                 bReadMolType = TRUE;
                             }
 
-                            push_molt(symtab, &nmol, molinfo, pline, wi);
+                            push_molt(symtab, molinfo, pline, wi);
+                            nmol = molinfo->size();
                             srenew(exclusionBlocks, nmol);
                             exclusionBlocks[nmol-1].nr      = 0;
-                            mi0                             = &((*molinfo)[nmol-1]);
+                            mi0                             = &molinfo->back();
                             mi0->atoms.haveMass             = TRUE;
                             mi0->atoms.haveCharge           = TRUE;
                             mi0->atoms.haveType             = TRUE;
@@ -791,7 +791,7 @@ static char **read_topol(const char *infile, const char *outfile,
                             int      whichmol;
                             bool     bCouple;
 
-                            push_mol(nmol, *molinfo, pline, &whichmol, &nrcopies, wi);
+                            push_mol(*molinfo, pline, &whichmol, &nrcopies, wi);
                             mi0 = &((*molinfo)[whichmol]);
                             molblock->resize(molblock->size() + 1);
                             molblock->back().type = whichmol;
@@ -924,7 +924,7 @@ static char **read_topol(const char *infile, const char *outfile,
     }
 
     DS_Done (&DS);
-    for (i = 0; i < nmol; i++)
+    for (int i = 0; i < gmx::index(molinfo->size()); i++)
     {
         gmx::doneExclusionBlocks(&(exclusionBlocks[i]));
     }
@@ -934,7 +934,7 @@ static char **read_topol(const char *infile, const char *outfile,
 
     if (*intermolecular_interactions != nullptr)
     {
-        sfree(mi0->atoms.atom);
+        sfree(intermolecular_interactions->get()->atoms.atom);
     }
 
     *nrmols = nmol;
@@ -942,24 +942,24 @@ static char **read_topol(const char *infile, const char *outfile,
     return title;
 }
 
-char **do_top(bool                          bVerbose,
-              const char                   *topfile,
-              const char                   *topppfile,
-              t_gromppopts                 *opts,
-              bool                          bZero,
-              t_symtab                     *symtab,
-              t_params                      plist[],
-              int                          *combination_rule,
-              double                       *repulsion_power,
-              real                         *fudgeQQ,
-              gpp_atomtype                 *atype,
-              int                          *nrmols,
-              t_molinfo                   **molinfo,
-              t_molinfo                   **intermolecular_interactions,
-              const t_inputrec             *ir,
-              std::vector<gmx_molblock_t>  *molblock,
-              bool                         *ffParametrizedWithHBondConstraints,
-              warninp                      *wi)
+char **do_top(bool                                  bVerbose,
+              const char                           *topfile,
+              const char                           *topppfile,
+              t_gromppopts                         *opts,
+              bool                                  bZero,
+              t_symtab                             *symtab,
+              t_params                              plist[],
+              int                                  *combination_rule,
+              double                               *repulsion_power,
+              real                                 *fudgeQQ,
+              gpp_atomtype                         *atype,
+              int                                  *nrmols,
+              std::vector<MoleculeInformation>     *molinfo,
+              std::unique_ptr<MoleculeInformation> *intermolecular_interactions,
+              const t_inputrec                     *ir,
+              std::vector<gmx_molblock_t>          *molblock,
+              bool                                 *ffParametrizedWithHBondConstraints,
+              warninp                              *wi)
 {
     /* Tmpfile might contain a long path */
     const char *tmpfile;
index c3ec5d5ff35e5c802a4180fbb70394f7fc49880c..c71d5b7a38e9187d5770b89d7e6252a680dfb590 100644 (file)
@@ -38,6 +38,7 @@
 #ifndef GMX_GMXPREPROCESS_TOPIO_H
 #define GMX_GMXPREPROCESS_TOPIO_H
 
+#include <memory>
 #include <vector>
 
 #include "gromacs/utility/real.h"
@@ -47,7 +48,7 @@ struct gmx_mtop_t;
 struct gpp_atomtype;
 struct t_gromppopts;
 struct t_inputrec;
-struct t_molinfo;
+struct MoleculeInformation;
 struct t_params;
 struct t_symtab;
 struct warninp;
@@ -57,24 +58,24 @@ typedef warninp *warninp_t;
 double check_mol(const gmx_mtop_t *mtop, warninp_t wi);
 /* Check mass and charge */
 
-char **do_top(bool                          bVerbose,
-              const char                   *topfile,
-              const char                   *topppfile,
-              t_gromppopts                 *opts,
-              bool                          bZero,
-              t_symtab                     *symtab,
-              t_params                      plist[],
-              int                          *combination_rule,
-              double                       *repulsion_power,
-              real                         *fudgeQQ,
-              gpp_atomtype                 *atype,
-              int                          *nrmols,
-              t_molinfo                   **molinfo,
-              t_molinfo                   **intermolecular_interactions,
-              const t_inputrec             *ir,
-              std::vector<gmx_molblock_t>  *molblock,
-              bool                         *ffParametrizedWithHBondConstraints,
-              warninp_t                     wi);
+char **do_top(bool                                               bVerbose,
+              const char                                        *topfile,
+              const char                                        *topppfile,
+              t_gromppopts                                      *opts,
+              bool                                               bZero,
+              t_symtab                                          *symtab,
+              t_params                                           plist[],
+              int                                               *combination_rule,
+              double                                            *repulsion_power,
+              real                                              *fudgeQQ,
+              gpp_atomtype                                      *atype,
+              int                                               *nrmols,
+              std::vector<MoleculeInformation>                  *molinfo,
+              std::unique_ptr<MoleculeInformation>              *intermolecular_interactions,
+              const t_inputrec                                  *ir,
+              std::vector<gmx_molblock_t>                       *molblock,
+              bool                                              *ffParametrizedWithHBondConstraints,
+              warninp_t                                          wi);
 
 /* This routine expects sys->molt[m].ilist to be of size F_NRE and ordered. */
 void generate_qmexcl(gmx_mtop_t *sys, t_inputrec *ir, warninp_t wi, GmxQmmmMode qmmmMode);
index abf5722ddfc16cdebb888604a5e9d1a3f7073f79..738a0afb9a0957d12e5adcf25f74af2042869d30 100644 (file)
@@ -1437,12 +1437,13 @@ void push_atom(t_symtab *symtab, t_block *cgs,
                   typeB == type ? ctype : ctypeB, mB, qB, wi);
 }
 
-void push_molt(t_symtab *symtab, int *nmol, t_molinfo **mol, char *line,
-               warninp *wi)
+void push_molt(t_symtab                         *symtab,
+               std::vector<MoleculeInformation> *mol,
+               char                             *line,
+               warninp                          *wi)
 {
     char       type[STRLEN];
-    int        nrexcl, i;
-    t_molinfo *newmol;
+    int        nrexcl;
 
     if ((sscanf(line, "%s%d", type, &nrexcl)) != 2)
     {
@@ -1450,26 +1451,22 @@ void push_molt(t_symtab *symtab, int *nmol, t_molinfo **mol, char *line,
     }
 
     /* Test if this moleculetype overwrites another */
-    i    = 0;
-    while (i < *nmol)
+    const auto found = std::find_if(mol->begin(), mol->end(),
+                                    [&type](const auto &m)
+                                    { return strcmp(*(m.name), type) == 0; });
+    if (found != mol->end())
     {
-        if (strcmp(*((*mol)[i].name), type) == 0)
-        {
-            auto message = gmx::formatString("moleculetype %s is redefined", type);
-            warning_error_and_exit(wi, message, FARGS);
-        }
-        i++;
+        auto message = gmx::formatString("moleculetype %s is redefined", type);
+        warning_error_and_exit(wi, message, FARGS);
     }
 
-    (*nmol)++;
-    srenew(*mol, *nmol);
-    newmol = &((*mol)[*nmol-1]);
-    init_molinfo(newmol);
+    mol->emplace_back();
+    mol->back().initMolInfo();
 
     /* Fill in the values */
-    newmol->name     = put_symtab(symtab, type);
-    newmol->nrexcl   = nrexcl;
-    newmol->excl_set = FALSE;
+    mol->back().name     = put_symtab(symtab, type);
+    mol->back().nrexcl   = nrexcl;
+    mol->back().excl_set = false;
 }
 
 static bool default_nb_params(int ftype, t_params bt[], t_atoms *at,
@@ -2371,7 +2368,7 @@ void push_vsitesn(Directive d, t_params bond[],
     sfree(weight);
 }
 
-void push_mol(int nrmols, t_molinfo mols[], char *pline, int *whichmol,
+void push_mol(gmx::ArrayRef<MoleculeInformation> mols, char *pline, int *whichmol,
               int *nrcopies,
               warninp *wi)
 {
@@ -2393,18 +2390,20 @@ void push_mol(int nrmols, t_molinfo mols[], char *pline, int *whichmol,
     int nrci    = 0;
     int matchci = -1;
     int matchcs = -1;
-    for (int i = 0; i < nrmols; i++)
+    int i       = 0;
+    for (const auto &mol : mols)
     {
-        if (strcmp(type, *(mols[i].name)) == 0)
+        if (strcmp(type, *(mol.name)) == 0)
         {
             nrcs++;
             matchcs = i;
         }
-        if (gmx_strcasecmp(type, *(mols[i].name)) == 0)
+        if (gmx_strcasecmp(type, *(mol.name)) == 0)
         {
             nrci++;
             matchci = i;
         }
+        i++;
     }
 
     if (nrcs == 1)
@@ -2576,7 +2575,7 @@ static void convert_pairs_to_pairsQ(t_params *plist,
     plist[F_LJ14].param = nullptr;
 }
 
-static void generate_LJCpairsNB(t_molinfo *mol, int nb_funct, t_params *nbp, warninp *wi)
+static void generate_LJCpairsNB(MoleculeInformation *mol, int nb_funct, t_params *nbp, warninp *wi)
 {
     int       n, ntype, i, j, k;
     t_atom   *atom;
@@ -2696,7 +2695,7 @@ static void decouple_atoms(t_atoms *atoms, int atomtype_decouple,
     }
 }
 
-void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple, real fudgeQQ,
+void convert_moltype_couple(MoleculeInformation *mol, int atomtype_decouple, real fudgeQQ,
                             int couple_lam0, int couple_lam1,
                             bool bCoupleIntra, int nb_funct, t_params *nbp,
                             warninp *wi)
index ee76fe3f00d2d44baff4f071c29aab0610a6940b..aa8d8b9bdb6789ecd746a0a5beef878972dd6e60 100644 (file)
@@ -38,6 +38,9 @@
 #ifndef GMX_GMXPREPROCESS_TOPPUSH_H
 #define GMX_GMXPREPROCESS_TOPPUSH_H
 
+#include <vector>
+
+#include "gromacs/utility/arrayref.h"
 #include "gromacs/utility/real.h"
 
 enum class Directive : int;
@@ -45,7 +48,7 @@ struct gpp_atomtype;
 struct gpp_bond_atomtype;
 struct t_atoms;
 struct t_block;
-struct t_molinfo;
+struct MoleculeInformation;
 struct t_nbparam;
 struct t_param;
 struct t_params;
@@ -104,11 +107,11 @@ void push_vsitesn(Directive d, t_params bond[],
                   t_atoms *at, char *line,
                   warninp *wi);
 
-void push_mol(int nrmols, t_molinfo mols[], char *pline,
+void push_mol(gmx::ArrayRef<MoleculeInformation> mols, char *pline,
               int *whichmol, int *nrcopies,
               warninp *wi);
 
-void push_molt(struct t_symtab *symtab, int *nmol, t_molinfo **mol, char *line,
+void push_molt(struct t_symtab *symtab, std::vector<MoleculeInformation> *mol, char *line,
                warninp *wi);
 
 void push_excl(char *line, gmx::ExclusionBlocks *b2, warninp *wi);
@@ -123,7 +126,7 @@ int add_atomtype_decoupled(struct t_symtab *symtab, gpp_atomtype *at,
  * Returns the atom type number.
  */
 
-void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple,
+void convert_moltype_couple(MoleculeInformation *mol, int atomtype_decouple,
                             real fudgeQQ,
                             int couple_lam0, int couple_lam1,
                             bool bCoupleIntra,
index c5a69f6014fd05e18cc241ab0a010617c4878c6a..ae1a78e340e50a8a177b82da2707c0d3eee4f857 100644 (file)
@@ -178,29 +178,6 @@ void add_param_to_list(t_params *list, t_param *b)
     list->nr++;
 }
 
-
-void init_molinfo(t_molinfo *mol)
-{
-    mol->nrexcl     = 0;
-    mol->excl_set   = FALSE;
-    mol->bProcessed = FALSE;
-    init_plist(mol->plist);
-    init_block(&mol->cgs);
-    init_block(&mol->mols);
-    init_blocka(&mol->excls);
-    init_t_atoms(&mol->atoms, 0, FALSE);
-}
-
-/* FREEING MEMORY */
-
-void done_mi(t_molinfo *mi)
-{
-    done_atom (&(mi->atoms));
-    done_block(&(mi->cgs));
-    done_block(&(mi->mols));
-    done_plist(mi->plist);
-}
-
 /* PRINTING STRUCTURES */
 
 static void print_bt(FILE *out, Directive d, gpp_atomtype *at,
index e581c812728fc8e9b6a0a96dd8558e3d17ded134..5664d606b11fe238d70dd95bef8ad5f639eab4f7 100644 (file)
@@ -45,7 +45,7 @@ struct gpp_atomtype;
 struct t_atoms;
 struct t_blocka;
 struct t_excls;
-struct t_molinfo;
+struct MoleculeInformation;
 struct t_param;
 struct t_params;
 
@@ -66,10 +66,6 @@ void add_param_to_list(t_params *list, t_param *b);
 void init_plist(t_params plist[]);
 void done_plist(t_params *plist);
 
-void init_molinfo(t_molinfo *mol);
-
-/* FREE */
-void done_mi(t_molinfo *mi);
 
 /* PRINTING */