Refactor preprocessing atom types.
[alexxy/gromacs.git] / src / gromacs / gmxpreprocess / gpp_atomtype.cpp
index c43d6bf0f9cd476b26bfc2679a571f60c8ca2b4e..8a478dd0c5e2f415a024dec8aeeeda86e990b8d6 100644 (file)
@@ -42,6 +42,8 @@
 #include <cmath>
 #include <cstring>
 
+#include <algorithm>
+
 #include "gromacs/gmxpreprocess/grompp_impl.h"
 #include "gromacs/gmxpreprocess/notset.h"
 #include "gromacs/gmxpreprocess/topdirs.h"
 #include "gromacs/utility/fatalerror.h"
 #include "gromacs/utility/smalloc.h"
 
-struct gpp_atomtype
+struct AtomTypeData
 {
-    int              nr;           /* The number of atomtypes          */
-    t_atom          *atom;         /* Array of atoms                   */
-    char          ***atomname;     /* Names of the atomtypes           */
-    t_param         *nb;           /* Nonbonded force default params   */
-    int             *bondatomtype; /* The bond_atomtype for each atomtype  */
-    int             *atomnumber;   /* Atomic number, used for QM/MM        */
+    //! Constructor initializes all fields.
+    AtomTypeData(const t_atom  &a,
+                 char         **name,
+                 const t_param *nb,
+                 const int      bondAtomType,
+                 const int      atomNumber) :
+        atom_(a), name_(name),
+        bondAtomType_(bondAtomType),
+        atomNumber_(atomNumber)
+    {
+        cp_param(&nb_, nb);
+    }
+    //! Actual atom data.
+    t_atom   atom_;
+    //! Atom name. The pointer is the result of a symtab operation and can be shallow copied.
+    char   **name_;
+    //! Nonbonded data.
+    t_param  nb_;
+    //! Bonded atomtype for the type.
+    int      bondAtomType_;
+    //! Atom number for the atom type.
+    int      atomNumber_;
 };
 
-int get_atomtype_type(const char *str, gpp_atomtype *ga)
+class PreprocessingAtomTypes::Impl
 {
-    int i;
+    public:
+        //! The number for currently loaded entries.
+        size_t size() const { return types.size(); }
+        //! The actual atom type data.
+        std::vector<AtomTypeData> types;
+};
 
+bool PreprocessingAtomTypes::isSet(int nt) const
+{
+    return ((nt >= 0) || (nt < gmx::ssize(*this)));
+}
+
+int PreprocessingAtomTypes::atomTypeFromName(const std::string &str) const
+{
     /* Atom types are always case sensitive */
-    for (i = 0; (i < ga->nr); i++)
+    auto found = std::find_if(impl_->types.begin(), impl_->types.end(),
+                              [&str](const auto &type)
+                              { return str == *type.name_; });
+    if (found == impl_->types.end())
     {
-        if (strcmp(str, *(ga->atomname[i])) == 0)
-        {
-            return i;
-        }
+        return NOTSET;
+    }
+    else
+    {
+        return std::distance(impl_->types.begin(), found);
     }
-
-    return NOTSET;
 }
 
-int get_atomtype_ntypes(gpp_atomtype *ga)
+size_t PreprocessingAtomTypes::size() const
 {
-    return ga->nr;
+    return impl_->size();
 }
 
-char *get_atomtype_name(int nt, gpp_atomtype *ga)
+const char *PreprocessingAtomTypes::atomNameFromAtomType(int nt) const
 {
-    if ((nt < 0) || (nt >= ga->nr))
-    {
-        return nullptr;
-    }
-
-    return *(ga->atomname[nt]);
+    return isSet(nt) ? *(impl_->types[nt].name_) : nullptr;
 }
 
-real get_atomtype_massA(int nt, gpp_atomtype *ga)
+real PreprocessingAtomTypes::atomMassAFromAtomType(int nt) const
 {
-    if ((nt < 0) || (nt >= ga->nr))
-    {
-        return NOTSET;
-    }
-
-    return ga->atom[nt].m;
+    return isSet(nt) ? impl_->types[nt].atom_.m : NOTSET;
 }
 
-real get_atomtype_massB(int nt, gpp_atomtype *ga)
+real PreprocessingAtomTypes::atomMassBFromAtomType(int nt) const
 {
-    if ((nt < 0) || (nt >= ga->nr))
-    {
-        return NOTSET;
-    }
-
-    return ga->atom[nt].mB;
+    return isSet(nt) ? impl_->types[nt].atom_.mB : NOTSET;
 }
 
-real get_atomtype_qA(int nt, gpp_atomtype *ga)
+real PreprocessingAtomTypes::atomChargeAFromAtomType(int nt) const
 {
-    if ((nt < 0) || (nt >= ga->nr))
-    {
-        return NOTSET;
-    }
-
-    return ga->atom[nt].q;
+    return isSet(nt) ? impl_->types[nt].atom_.q : NOTSET;
 }
 
-real get_atomtype_qB(int nt, gpp_atomtype *ga)
+real PreprocessingAtomTypes::atomChargeBFromAtomType(int nt) const
 {
-    if ((nt < 0) || (nt >= ga->nr))
-    {
-        return NOTSET;
-    }
-
-    return ga->atom[nt].qB;
+    return isSet(nt) ? impl_->types[nt].atom_.qB : NOTSET;
 }
 
-int get_atomtype_ptype(int nt, gpp_atomtype *ga)
+int PreprocessingAtomTypes::atomParticleTypeFromAtomType(int nt) const
 {
-    if ((nt < 0) || (nt >= ga->nr))
-    {
-        return NOTSET;
-    }
-
-    return ga->atom[nt].ptype;
+    return isSet(nt) ? impl_->types[nt].atom_.ptype : NOTSET;
 }
 
-int get_atomtype_batype(int nt, const gpp_atomtype* ga)
+int PreprocessingAtomTypes::bondAtomTypeFromAtomType(int nt) const
 {
-    if ((nt < 0) || (nt >= ga->nr))
-    {
-        return NOTSET;
-    }
-
-    return ga->bondatomtype[nt];
+    return isSet(nt) ? impl_->types[nt].bondAtomType_ : NOTSET;
 }
 
-int get_atomtype_atomnumber(int nt, gpp_atomtype *ga)
+int PreprocessingAtomTypes::atomNumberFromAtomType(int nt) const
 {
-    if ((nt < 0) || (nt >= ga->nr))
-    {
-        return NOTSET;
-    }
-
-    return ga->atomnumber[nt];
+    return isSet(nt) ? impl_->types[nt].atomNumber_ : NOTSET;
 }
 
-real get_atomtype_nbparam(int nt, int param, gpp_atomtype *ga)
+real PreprocessingAtomTypes::atomNonBondedParamFromAtomType(int nt, int param) const
 {
-    if ((nt < 0) || (nt >= ga->nr))
+    if (!isSet(nt))
     {
         return NOTSET;
     }
@@ -174,113 +166,101 @@ real get_atomtype_nbparam(int nt, int param, gpp_atomtype *ga)
     {
         return NOTSET;
     }
-    return ga->nb[nt].c[param];
+    return impl_->types[nt].nb_.c[param];
 }
 
-gpp_atomtype *init_atomtype()
-{
-    gpp_atomtype *ga;
-
-    snew(ga, 1);
+PreprocessingAtomTypes::PreprocessingAtomTypes()
+    : impl_(new Impl)
+{}
 
-    ga->nr           = 0;
-    ga->atom         = nullptr;
-    ga->atomname     = nullptr;
-    ga->nb           = nullptr;
-    ga->bondatomtype = nullptr;
-    ga->atomnumber   = nullptr;
-
-    return ga;
-}
+PreprocessingAtomTypes::PreprocessingAtomTypes(PreprocessingAtomTypes &&old) noexcept
+    : impl_(std::move(old.impl_))
+{}
 
-int set_atomtype(int nt, gpp_atomtype *ga, t_symtab *tab,
-                 t_atom *a, const char *name, t_param *nb,
-                 int bondatomtype, int atomnumber)
+PreprocessingAtomTypes &PreprocessingAtomTypes::operator=(PreprocessingAtomTypes &&old) noexcept
 {
-    if ((nt < 0) || (nt >= ga->nr))
-    {
-        return NOTSET;
-    }
-
-    ga->atom[nt]         = *a;
-    ga->atomname[nt]     = put_symtab(tab, name);
-    ga->nb[nt]           = *nb;
-    ga->bondatomtype[nt] = bondatomtype;
-    ga->atomnumber[nt]   = atomnumber;
-
-    return nt;
+    impl_ = std::move(old.impl_);
+    return *this;
 }
 
-int add_atomtype(gpp_atomtype *ga, t_symtab *tab,
-                 t_atom *a, const char *name, t_param *nb,
-                 int bondatomtype, int atomnumber)
+PreprocessingAtomTypes::~PreprocessingAtomTypes()
+{}
+
+int PreprocessingAtomTypes::addType(t_symtab          *tab,
+                                    const t_atom      &a,
+                                    const char        *name,
+                                    const t_param     *nb,
+                                    int                bondAtomType,
+                                    int                atomNumber)
 {
-    int i;
+    auto found = std::find_if(impl_->types.begin(), impl_->types.end(),
+                              [&name](const AtomTypeData &data)
+                              { return strcmp(name, *data.name_) == 0; });
 
-    for (i = 0; (i < ga->nr); i++)
-    {
-        if (strcmp(*ga->atomname[i], name) == 0)
-        {
-            break;
-        }
-    }
-    if (i == ga->nr)
+    if (found == impl_->types.end())
     {
-        ga->nr++;
-        srenew(ga->atom, ga->nr);
-        srenew(ga->atomname, ga->nr);
-        srenew(ga->nb, ga->nr);
-        srenew(ga->bondatomtype, ga->nr);
-        srenew(ga->atomnumber, ga->nr);
-
-        return set_atomtype(ga->nr-1, ga, tab, a, name, nb, bondatomtype, atomnumber);
+        impl_->types.emplace_back(a,
+                                  put_symtab(tab, name),
+                                  nb,
+                                  bondAtomType,
+                                  atomNumber);
+        return size() - 1;
     }
     else
     {
-        return i;
+        return std::distance(impl_->types.begin(), found);
     }
 }
 
-void print_at (FILE * out, gpp_atomtype *ga)
+int PreprocessingAtomTypes::setType(int                nt,
+                                    t_symtab          *tab,
+                                    const t_atom      &a,
+                                    const char        *name,
+                                    const t_param     *nb,
+                                    int                bondAtomType,
+                                    int                atomNumber)
 {
-    int         i;
-    t_atom     *atom = ga->atom;
-    t_param    *nb   = ga->nb;
+    if (!isSet(nt))
+    {
+        return NOTSET;
+    }
+
+    impl_->types[nt].atom_         = a;
+    impl_->types[nt].name_         = put_symtab(tab, name);
+    cp_param(&impl_->types[nt].nb_, nb);
+    impl_->types[nt].bondAtomType_ = bondAtomType;
+    impl_->types[nt].atomNumber_   = atomNumber;
+
+    return nt;
+}
 
+void PreprocessingAtomTypes::printTypes(FILE * out)
+{
     fprintf (out, "[ %s ]\n", dir2str(Directive::d_atomtypes));
     fprintf (out, "; %6s  %8s  %8s  %8s  %12s  %12s\n",
              "type", "mass", "charge", "particle", "c6", "c12");
-    for (i = 0; (i < ga->nr); i++)
+    for (auto &entry : impl_->types)
     {
         fprintf(out, "%8s  %8.3f  %8.3f  %8s  %12e  %12e\n",
-                *(ga->atomname[i]), atom[i].m, atom[i].q, "A",
-                nb[i].c0(), nb[i].c1());
+                *(entry.name_), entry.atom_.m, entry.atom_.q, "A",
+                entry.nb_.c0(), entry.nb_.c1());
     }
 
     fprintf (out, "\n");
 }
 
-void done_atomtype(gpp_atomtype *ga)
+static int search_atomtypes(const PreprocessingAtomTypes *ga,
+                            int                          *n,
+                            gmx::ArrayRef<int>            typelist,
+                            int                           thistype,
+                            t_param                       param[],
+                            int                           ftype)
 {
-    sfree(ga->atom);
-    sfree(ga->atomname);
-    sfree(ga->nb);
-    sfree(ga->bondatomtype);
-    sfree(ga->atomnumber);
-    ga->nr = 0;
-    sfree(ga);
-}
-
-static int search_atomtypes(gpp_atomtype *ga, int *n, int typelist[],
-                            int thistype,
-                            t_param param[], int ftype)
-{
-    int      i, nn, nrfp, j, k, ntype, tli;
-    bool     bFound = FALSE;
+    int      i, nn, nrfp, ntype;
 
     nn    = *n;
     nrfp  = NRFP(ftype);
-    ntype = get_atomtype_ntypes(ga);
+    ntype = ga->size();
 
     for (i = 0; (i < nn); i++)
     {
@@ -292,19 +272,19 @@ static int search_atomtypes(gpp_atomtype *ga, int *n, int typelist[],
 
         /* Otherwise, check if the parameters are identical to any previously added type */
 
-        bFound = TRUE;
-        for (j = 0; j < ntype && bFound; j++)
+        bool bFound = true;
+        for (int j = 0; j < ntype && bFound; j++)
         {
             /* Check nonbonded parameters */
-            for (k = 0; k < nrfp && bFound; k++)
+            for (int k = 0; k < nrfp && bFound; k++)
             {
                 bFound = (param[ntype*typelist[i]+j].c[k] == param[ntype*thistype+j].c[k]);
             }
 
             /* Check atomnumber */
-            tli    = typelist[i];
+            int tli    = typelist[i];
             bFound = bFound &&
-                (get_atomtype_atomnumber(tli, ga) == get_atomtype_atomnumber(thistype, ga));
+                (ga->atomNumberFromAtomType(tli) == ga->atomNumberFromAtomType(thistype));
         }
         if (bFound)
         {
@@ -326,19 +306,15 @@ static int search_atomtypes(gpp_atomtype *ga, int *n, int typelist[],
     return i;
 }
 
-void renum_atype(gmx::ArrayRef<InteractionTypeParameters> plist, gmx_mtop_t *mtop,
-                 int *wall_atomtype,
-                 gpp_atomtype *ga, bool bVerbose)
+void PreprocessingAtomTypes::renumberTypes(gmx::ArrayRef<InteractionTypeParameters> plist,
+                                           gmx_mtop_t                              *mtop,
+                                           int                                     *wall_atomtype,
+                                           bool                                     bVerbose)
 {
-    int         i, j, k, l, mi, mj, nat, nrfp, ftype, ntype;
-    t_atoms    *atoms;
-    t_param    *nbsnew;
-    int        *typelist;
-    int        *new_atomnumber;
-    char     ***new_atomname;
+    int         nat, ftype, ntype;
 
-    ntype = get_atomtype_ntypes(ga);
-    snew(typelist, ntype);
+    ntype = size();
+    std::vector<int> typelist(ntype);
 
     if (bVerbose)
     {
@@ -368,57 +344,58 @@ void renum_atype(gmx::ArrayRef<InteractionTypeParameters> plist, gmx_mtop_t *mto
      * can determine if two types should be merged.
      */
     nat = 0;
-    for (gmx_moltype_t &moltype : mtop->moltype)
+    for (const gmx_moltype_t &moltype : mtop->moltype)
     {
-        atoms = &moltype.atoms;
-        for (i = 0; (i < atoms->nr); i++)
+        const t_atoms *atoms = &moltype.atoms;
+        for (int i = 0; (i < atoms->nr); i++)
         {
             atoms->atom[i].type =
-                search_atomtypes(ga, &nat, typelist, atoms->atom[i].type,
+                search_atomtypes(this, &nat, typelist, atoms->atom[i].type,
                                  plist[ftype].param, ftype);
             atoms->atom[i].typeB =
-                search_atomtypes(ga, &nat, typelist, atoms->atom[i].typeB,
+                search_atomtypes(this, &nat, typelist, atoms->atom[i].typeB,
                                  plist[ftype].param, ftype);
         }
     }
 
-    for (i = 0; i < 2; i++)
+    for (int i = 0; i < 2; i++)
     {
         if (wall_atomtype[i] >= 0)
         {
-            wall_atomtype[i] = search_atomtypes(ga, &nat, typelist, wall_atomtype[i],
+            wall_atomtype[i] = search_atomtypes(this, &nat, typelist, wall_atomtype[i],
                                                 plist[ftype].param, ftype);
         }
     }
 
-    snew(new_atomnumber, nat);
-    snew(new_atomname, nat);
+    std::vector<AtomTypeData> new_types;
     /* We now have a list of unique atomtypes in typelist */
 
     /* Renumber nlist */
-    nbsnew = nullptr;
+    /* Renumber nlist */
+    t_param *nbsnew = nullptr;
     snew(nbsnew, plist[ftype].nr);
 
-    nrfp  = NRFP(ftype);
+    int nrfp  = NRFP(ftype);
 
-    for (i = k = 0; (i < nat); i++)
+    int k = 0;
+    for (int i = 0; (i < nat); i++)
     {
-        mi = typelist[i];
-        for (j = 0; (j < nat); j++, k++)
+        int mi = typelist[i];
+        for (int j = 0; (j < nat); j++, k++)
         {
-            mj = typelist[j];
-            for (l = 0; (l < nrfp); l++)
+            int mj = typelist[j];
+            for (int l = 0; (l < nrfp); l++)
             {
                 nbsnew[k].c[l] = plist[ftype].param[ntype*mi+mj].c[l];
             }
         }
-        new_atomnumber[i] = get_atomtype_atomnumber(mi, ga);
-        new_atomname[i]   = ga->atomname[mi];
+        new_types.push_back(impl_->types[mi]);
     }
 
+    int i;
     for (i = 0; (i < nat*nat); i++)
     {
-        for (l = 0; (l < nrfp); l++)
+        for (int l = 0; (l < nrfp); l++)
         {
             plist[ftype].param[i].c[l] = nbsnew[i].c[l];
         }
@@ -426,30 +403,20 @@ void renum_atype(gmx::ArrayRef<InteractionTypeParameters> plist, gmx_mtop_t *mto
     plist[ftype].nr     = i;
     mtop->ffparams.atnr = nat;
 
-    sfree(ga->atomnumber);
-    /* Dangling atomname pointers ? */
-    sfree(ga->atomname);
-
-    ga->atomnumber = new_atomnumber;
-    ga->atomname   = new_atomname;
-
-    ga->nr = nat;
+    impl_->types = new_types;
 
     sfree(nbsnew);
-    sfree(typelist);
 }
 
-void copy_atomtype_atomtypes(gpp_atomtype *ga, t_atomtypes *atomtypes)
+void PreprocessingAtomTypes::copyTot_atomtypes(t_atomtypes *atomtypes) const
 {
-    int i, ntype;
-
     /* Copy the atomtype data to the topology atomtype list */
-    ntype         = get_atomtype_ntypes(ga);
+    int ntype         = size();
     atomtypes->nr = ntype;
     snew(atomtypes->atomnumber, ntype);
 
-    for (i = 0; i < ntype; i++)
+    for (int i = 0; i < ntype; i++)
     {
-        atomtypes->atomnumber[i] = ga->atomnumber[i];
+        atomtypes->atomnumber[i] = impl_->types[i].atomNumber_;
     }
 }