From 2fadcf2debf5bf2eecfa3911c4959e3f2ea06619 Mon Sep 17 00:00:00 2001 From: Paul Bauer Date: Fri, 25 Jan 2019 17:07:02 +0100 Subject: [PATCH] Refactor preprocessing atom types. Refs #2833 Change-Id: Ifd7f583fd5908c1ce225e379b58757f9a09b2100 --- src/gromacs/gmxpreprocess/convparm.cpp | 11 +- src/gromacs/gmxpreprocess/gen_vsite.cpp | 14 +- src/gromacs/gmxpreprocess/gen_vsite.h | 4 +- src/gromacs/gmxpreprocess/gpp_atomtype.cpp | 369 +++++++++++++---------------- src/gromacs/gmxpreprocess/gpp_atomtype.h | 231 ++++++++++++++---- src/gromacs/gmxpreprocess/grompp.cpp | 38 ++- src/gromacs/gmxpreprocess/nm2type.cpp | 62 ++--- src/gromacs/gmxpreprocess/nm2type.h | 4 +- src/gromacs/gmxpreprocess/pdb2gmx.cpp | 11 +- src/gromacs/gmxpreprocess/pdb2top.cpp | 4 +- src/gromacs/gmxpreprocess/pdb2top.h | 6 +- src/gromacs/gmxpreprocess/resall.cpp | 23 +- src/gromacs/gmxpreprocess/resall.h | 18 +- src/gromacs/gmxpreprocess/ter_db.cpp | 22 +- src/gromacs/gmxpreprocess/ter_db.h | 4 +- src/gromacs/gmxpreprocess/tomorse.cpp | 10 +- src/gromacs/gmxpreprocess/tomorse.h | 4 +- src/gromacs/gmxpreprocess/topio.cpp | 97 ++++---- src/gromacs/gmxpreprocess/topio.h | 37 ++- src/gromacs/gmxpreprocess/toppush.cpp | 190 ++++++++------- src/gromacs/gmxpreprocess/toppush.h | 32 +-- src/gromacs/gmxpreprocess/topshake.cpp | 42 ++-- src/gromacs/gmxpreprocess/toputil.cpp | 69 +++--- src/gromacs/gmxpreprocess/toputil.h | 8 +- src/gromacs/gmxpreprocess/vsite_parm.cpp | 34 ++- src/gromacs/gmxpreprocess/vsite_parm.h | 4 +- src/gromacs/gmxpreprocess/x2top.cpp | 47 ++-- 27 files changed, 739 insertions(+), 656 deletions(-) diff --git a/src/gromacs/gmxpreprocess/convparm.cpp b/src/gromacs/gmxpreprocess/convparm.cpp index 63f7aca8f8..546cd97ae0 100644 --- a/src/gromacs/gmxpreprocess/convparm.cpp +++ b/src/gromacs/gmxpreprocess/convparm.cpp @@ -109,11 +109,10 @@ static int assign_param(t_functype ftype, t_iparams *newparam, real old[MAXFORCEPARAM], int comb, double reppow) { - int i, j; - bool all_param_zero = TRUE; + bool all_param_zero = true; /* Set to zero */ - for (j = 0; (j < MAXFORCEPARAM); j++) + for (int j = 0; (j < MAXFORCEPARAM); j++) { newparam->generic.buf[j] = 0.0; /* If all parameters are zero we might not add some interaction types (selected below). @@ -194,7 +193,7 @@ assign_param(t_functype ftype, t_iparams *newparam, break; case F_QUARTIC_ANGLES: newparam->qangle.theta = old[0]; - for (i = 0; i < 5; i++) + for (int i = 0; i < 5; i++) { newparam->qangle.c[i] = old[i+1]; } @@ -365,14 +364,14 @@ assign_param(t_functype ftype, t_iparams *newparam, newparam->dihres.kfacB = old[5]; break; case F_RBDIHS: - for (i = 0; (i < NR_RBDIHS); i++) + for (int i = 0; (i < NR_RBDIHS); i++) { newparam->rbdihs.rbcA[i] = old[i]; newparam->rbdihs.rbcB[i] = old[NR_RBDIHS+i]; } break; case F_CBTDIHS: - for (i = 0; (i < NR_CBTDIHS); i++) + for (int i = 0; (i < NR_CBTDIHS); i++) { newparam->cbtdihs.cbtcA[i] = old[i]; } diff --git a/src/gromacs/gmxpreprocess/gen_vsite.cpp b/src/gromacs/gmxpreprocess/gen_vsite.cpp index c5ce5f54db..1a50d6b100 100644 --- a/src/gromacs/gmxpreprocess/gen_vsite.cpp +++ b/src/gromacs/gmxpreprocess/gen_vsite.cpp @@ -570,11 +570,11 @@ static int get_atype(int atom, t_atoms *at, gmx::ArrayRefatomTypeFromName(name); if (tp == NOTSET) { gmx_fatal(FARGS, "Dummy mass type (%s) not found in atom type database", @@ -871,7 +871,7 @@ static void calc_vsite3_param(real xd, real yd, real xi, real yi, real xj, real } -static int gen_vsites_trp(gpp_atomtype *atype, +static int gen_vsites_trp(PreprocessingAtomTypes *atype, std::vector *newx, t_atom *newatom[], char ***newatomname[], int *o2n[], int *newvsite_type[], int *newcgnr[], @@ -1145,7 +1145,7 @@ static int gen_vsites_trp(gpp_atomtype *atype, } -static int gen_vsites_tyr(gpp_atomtype *atype, +static int gen_vsites_tyr(PreprocessingAtomTypes *atype, std::vector *newx, t_atom *newatom[], char ***newatomname[], int *o2n[], int *newvsite_type[], int *newcgnr[], @@ -1549,7 +1549,7 @@ static bool is_vsite(int vsite_type) static char atomnamesuffix[] = "1234"; -void do_vsites(gmx::ArrayRef rtpFFDB, gpp_atomtype *atype, +void do_vsites(gmx::ArrayRef rtpFFDB, PreprocessingAtomTypes *atype, t_atoms *at, t_symtab *symtab, std::vector *x, gmx::ArrayRef plist, int *vsite_type[], int *cgnr[], @@ -1790,7 +1790,7 @@ void do_vsites(gmx::ArrayRef rtpFFDB, gpp_atomtype *aty &nrbonds, &nrHatoms, Hatoms, &Heavy, &nrheavies, heavies); /* get Heavy atom type */ tpHeavy = get_atype(Heavy, at, rtpFFDB, &rt); - strcpy(tpname, get_atomtype_name(tpHeavy, atype)); + strcpy(tpname, atype->atomNameFromAtomType(tpHeavy)); bWARNING = FALSE; bAddVsiteParam = TRUE; @@ -1896,7 +1896,7 @@ void do_vsites(gmx::ArrayRef rtpFFDB, gpp_atomtype *aty } /* get dummy mass type from first char of heavy atom type (N or C) */ - strcpy(nexttpname, get_atomtype_name(get_atype(heavies[0], at, rtpFFDB, &rt), atype)); + strcpy(nexttpname, atype->atomNameFromAtomType(get_atype(heavies[0], at, rtpFFDB, &rt))); std::string ch = get_dummymass_name(vsiteconflist, tpname, nexttpname); std::string name; if (ch.empty()) diff --git a/src/gromacs/gmxpreprocess/gen_vsite.h b/src/gromacs/gmxpreprocess/gen_vsite.h index a5d89a1f11..532d6b79e6 100644 --- a/src/gromacs/gmxpreprocess/gen_vsite.h +++ b/src/gromacs/gmxpreprocess/gen_vsite.h @@ -42,7 +42,7 @@ #include "gromacs/utility/arrayref.h" #include "gromacs/utility/real.h" -struct gpp_atomtype; +class PreprocessingAtomTypes; struct t_atoms; struct InteractionTypeParameters; struct PreprocessResidue; @@ -50,7 +50,7 @@ struct t_symtab; /* stuff for pdb2gmx */ -void do_vsites(gmx::ArrayRef rtpFFDB, gpp_atomtype *atype, +void do_vsites(gmx::ArrayRef rtpFFDB, PreprocessingAtomTypes *atype, t_atoms *at, t_symtab *symtab, std::vector *x, gmx::ArrayRef plist, int *dummy_type[], int *cgnr[], real mHmult, bool bVSiteAromatics, diff --git a/src/gromacs/gmxpreprocess/gpp_atomtype.cpp b/src/gromacs/gmxpreprocess/gpp_atomtype.cpp index c43d6bf0f9..8a478dd0c5 100644 --- a/src/gromacs/gmxpreprocess/gpp_atomtype.cpp +++ b/src/gromacs/gmxpreprocess/gpp_atomtype.cpp @@ -42,6 +42,8 @@ #include #include +#include + #include "gromacs/gmxpreprocess/grompp_impl.h" #include "gromacs/gmxpreprocess/notset.h" #include "gromacs/gmxpreprocess/topdirs.h" @@ -53,120 +55,110 @@ #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 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 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 plist, gmx_mtop_t *mtop, - int *wall_atomtype, - gpp_atomtype *ga, bool bVerbose) +void PreprocessingAtomTypes::renumberTypes(gmx::ArrayRef 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 typelist(ntype); if (bVerbose) { @@ -368,57 +344,58 @@ void renum_atype(gmx::ArrayRef 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 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 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_; } } diff --git a/src/gromacs/gmxpreprocess/gpp_atomtype.h b/src/gromacs/gmxpreprocess/gpp_atomtype.h index 075f7c895a..57e0c95548 100644 --- a/src/gromacs/gmxpreprocess/gpp_atomtype.h +++ b/src/gromacs/gmxpreprocess/gpp_atomtype.h @@ -41,65 +41,194 @@ #include #include "gromacs/utility/arrayref.h" +#include "gromacs/utility/classhelpers.h" #include "gromacs/utility/real.h" struct gmx_mtop_t; -struct gpp_atomtype; struct t_atom; struct t_atomtypes; struct t_param; struct InteractionTypeParameters; struct t_symtab; -int get_atomtype_type(const char *str, gpp_atomtype *at); -/* Return atomtype corresponding to case-insensitive str - or NOTSET if not found */ - -int get_atomtype_ntypes(gpp_atomtype *at); -/* Return number of atomtypes */ - -char *get_atomtype_name(int nt, gpp_atomtype *at); -/* Return name corresponding to atomtype nt, or NULL if not found */ - -real get_atomtype_massA(int nt, gpp_atomtype *at); -real get_atomtype_massB(int nt, gpp_atomtype *at); -real get_atomtype_qA(int nt, gpp_atomtype *at); -real get_atomtype_qB(int nt, gpp_atomtype *at); -int get_atomtype_ptype(int nt, gpp_atomtype *at); -int get_atomtype_batype(int nt, const gpp_atomtype* at); -int get_atomtype_atomnumber(int nt, gpp_atomtype *at); - -/* Return the above variable for atomtype nt, or NOTSET if not found */ - -real get_atomtype_nbparam(int nt, int param, gpp_atomtype *at); -/* Similar to the previous but returns the paramth parameter or NOTSET */ - -gpp_atomtype *init_atomtype(); -/* Return a new atomtype structure */ - -void done_atomtype(gpp_atomtype *at); -/* Free the memory in the structure */ - -int set_atomtype(int nt, gpp_atomtype *at, t_symtab *tab, - t_atom *a, const char *name, t_param *nb, - int bondatomtype, int atomnumber); -/* Set the values of an existing atom type nt. Returns nt on success or - NOTSET on error. */ - -int add_atomtype(gpp_atomtype *at, t_symtab *tab, - t_atom *a, const char *name, t_param *nb, - int bondatomtype, int atomnumber); -/* Add a complete new atom type to an existing atomtype structure. Returns - the number of the atom type. */ - -void print_at (FILE * out, gpp_atomtype *at); -/* Print an atomtype record to a text file */ - -void renum_atype(gmx::ArrayRef plist, gmx_mtop_t *mtop, - int *wall_atomtype, - gpp_atomtype *at, bool bVerbose); - -void copy_atomtype_atomtypes(gpp_atomtype *atype, t_atomtypes *atypes); -/* Copy from one structure to another */ +/*! \libinternal \brief + * Storage of all atom types used during preprocessing of a simulation + * input. + */ +class PreprocessingAtomTypes +{ + public: + PreprocessingAtomTypes(); + //! Move constructor. + PreprocessingAtomTypes(PreprocessingAtomTypes &&old) noexcept; + //! Move assignment constructor. + PreprocessingAtomTypes &operator=(PreprocessingAtomTypes &&old) noexcept; + + ~PreprocessingAtomTypes(); + + /*! \brief + * Get atom type index for atom type name if present in the database, or NOTSET. + * + * \todo The code should be changed to instead use a gmx::compat version + * of std::optional to return an iterator to the element being searched, + * or an empty optional construct if the entry has not been found. + * + * \param[in] str Input string to search type for. + * \returns Atomtype as integer. + */ + int atomTypeFromName(const std::string &str) const; + + //! Get number of defined atom types. + size_t size() const; + + /*! \brief + * Get name of atom from internal atom type number. + * + * \param[in] nt Internal number of atom type. + * \returns The type name. + */ + const char *atomNameFromAtomType(int nt) const; + + /*! \brief + * Get normal mass of atom from internal atom type number. + * + * \param[in] nt Internal number of atom type. + * \returns The mass for the atom in state A or NOTSET. + */ + real atomMassAFromAtomType(int nt) const; + + /*! \brief + * Get mass for B state of atom from internal atom type number. + * + * \param[in] nt Internal number of atom type. + * \returns The mass for the atom in state B or NOTSET. + */ + real atomMassBFromAtomType(int nt) const; + + /*! \brief + * Get normal charge of atom from internal atom type number. + * + * \param[in] nt Internal number of atom type. + * \returns The charge for the atom in state A or NOTSET. + */ + real atomChargeAFromAtomType(int nt) const; + + /*! \brief + * Get charge for B state of atom from internal atom type number. + * + * \param[in] nt Internal number of atom type. + * \returns The charge for the atom in state B or NOTSET. + */ + real atomChargeBFromAtomType(int nt) const; + + /*! \brief + * Get particle type for atom type \p nt + * + * \param[in] nt Internal number of atom type. + * \returns The particle type or NOTSET. + */ + int atomParticleTypeFromAtomType(int nt) const; + + /*! \brief + * Get bond atom parameter of atom from internal atom type number. + * + * \param[in] nt Internal number of atom type. + * \returns The bond atom parameter or NOTSET. + */ + int bondAtomTypeFromAtomType(int nt) const; + + /*! \brief + * Get atomic number of atom from internal atom type number. + * + * \param[in] nt Internal number of atom type. + * \returns The atomic number type or NOTSET. + */ + int atomNumberFromAtomType(int nt) const; + + /*! \brief + * Get the value of \p param of type \p nt. + * + * \param[in] param The parameter value to find. + * \param[in] nt The number of the type. + * \returns The value of the parameter or NOTSET. + */ + real atomNonBondedParamFromAtomType(int nt, int param) const; + + /*! \brief + * If a value is within the range of the current types or not. + * + * \param[in] nt Value to check. + * \returns True if value is in range. + */ + bool isSet(int nt) const; + + /*! \brief + * Print data to file. + * + * \param[in] out File pointer. + */ + void printTypes(FILE *out); + + /*! \brief + * Set the values of an existing atom type \p nt. + * + * \param[in] nt Type that should be set. + * \param[in] tab Symbol table. + * \param[in] a Atom information. + * \param[in] name Atom name. + * \param[in] nb Nonbonded parameters. + * \param[in] bondAtomType What kind of bonded interactions are there. + * \param[in] atomNumber Atomic number of the entry. + * \returns Number of the type set or NOTSET + */ + int setType(int nt, + t_symtab *tab, + const t_atom &a, + const char *name, + const t_param *nb, + int bondAtomType, + int atomNumber); + + /*! \brief + * Add new atom type to database. + * + * \param[in] tab Symbol table. + * \param[in] a Atom information. + * \param[in] name Atom name. + * \param[in] nb Nonbonded parameters. + * \param[in] bondAtomType What kind of bonded interactions are there. + * \param[in] atomNumber Atomic number of the entry. + * \returns Number of entries in database. + */ + int addType(t_symtab *tab, + const t_atom &a, + const char *name, + const t_param *nb, + int bondAtomType, + int atomNumber); + + /*! \brief + * Renumber existing atom type entries. + * + * \param[in] plist List of parameters. + * \param[in] mtop Global topology. + * \param[inout] wallAtomType Atom types of wall atoms, which may also be renumbered + * \param[in] verbose If we want to print additional info. + */ + void renumberTypes(gmx::ArrayRef plist, + gmx_mtop_t *mtop, + int *wallAtomType, + bool verbose); + + /*! \brief + * Copy information to other structure. + * + * \param[inout] atypes Other datastructure to copy to. + */ + void copyTot_atomtypes(t_atomtypes *atypes) const; + private: + class Impl; + gmx::PrivateImplPointer impl_; +}; #endif diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index 24fadfbb0c..1d6430be8b 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -507,7 +507,7 @@ static void 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, + PreprocessingAtomTypes *atypes, gmx_mtop_t *sys, std::vector *mi, std::unique_ptr *intermolecular_interactions, gmx::ArrayRef plist, @@ -516,7 +516,7 @@ new_status(const char *topfile, const char *topppfile, const char *confin, warninp *wi) { std::vector molblock; - int i, nrmols, nmismatch; + int i, nmismatch; bool ffParametrizedWithHBondConstraints; char buf[STRLEN]; char warn_buf[STRLEN]; @@ -524,7 +524,7 @@ new_status(const char *topfile, const char *topppfile, const char *confin, /* TOPOLOGY processing */ sys->name = do_top(bVerbose, topfile, topppfile, opts, bZero, &(sys->symtab), plist, comb, reppow, fudgeQQ, - atype, &nrmols, mi, intermolecular_interactions, + atypes, mi, intermolecular_interactions, ir, &molblock, &ffParametrizedWithHBondConstraints, @@ -557,7 +557,7 @@ new_status(const char *topfile, const char *topppfile, const char *confin, if (bMorse) { - convert_harmonics(*mi, atype); + convert_harmonics(*mi, atypes); } if (ir->eDisre == edrNone) @@ -1012,7 +1012,7 @@ static void gen_posres(gmx_mtop_t *mtop, read_posres(mtop, mi, TRUE, fnB, rc_scaling, ePBC, comB, wi); } -static void set_wall_atomtype(gpp_atomtype *at, t_gromppopts *opts, +static void set_wall_atomtype(PreprocessingAtomTypes *at, t_gromppopts *opts, t_inputrec *ir, warninp *wi) { int i; @@ -1024,7 +1024,7 @@ static void set_wall_atomtype(gpp_atomtype *at, t_gromppopts *opts, } for (i = 0; i < ir->nwall; i++) { - ir->wall_atomtype[i] = get_atomtype_type(opts->wall_atomtype[i], at); + ir->wall_atomtype[i] = at->atomTypeFromName(opts->wall_atomtype[i]); if (ir->wall_atomtype[i] == NOTSET) { sprintf(warn_buf, "Specified wall atom type %s is not defined", opts->wall_atomtype[i]); @@ -1692,12 +1692,10 @@ int gmx_grompp(int argc, char *argv[]) t_gromppopts *opts; std::vector mi; std::unique_ptr intermolecular_interactions; - gpp_atomtype *atype; int nvsite, comb; real fudgeQQ; double reppow; const char *mdparin; - int ntype; bool bNeedVel, bGenVel; gmx_bool have_atomnumber; gmx_output_env_t *oenv; @@ -1800,8 +1798,8 @@ int gmx_grompp(int argc, char *argv[]) std::array plist; init_plist(plist); - gmx_mtop_t sys; - atype = init_atomtype(); + gmx_mtop_t sys; + PreprocessingAtomTypes atypes; if (debug) { pr_symtab(debug, 0, "Just opened", &sys.symtab); @@ -1816,7 +1814,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, &mi, &intermolecular_interactions, + &atypes, &sys, &mi, &intermolecular_interactions, plist, &comb, &reppow, &fudgeQQ, opts->bMorse, wi); @@ -1831,7 +1829,7 @@ int gmx_grompp(int argc, char *argv[]) for (size_t mt = 0; mt < sys.moltype.size(); mt++) { nvsite += - set_vsites(bVerbose, &sys.moltype[mt].atoms, atype, mi[mt].plist); + set_vsites(bVerbose, &sys.moltype[mt].atoms, &atypes, mi[mt].plist); } /* now throw away all obsolete bonds, angles and dihedrals: */ /* note: constraints are ALWAYS removed */ @@ -1875,9 +1873,9 @@ int gmx_grompp(int argc, char *argv[]) /* If we are doing QM/MM, check that we got the atom numbers */ have_atomnumber = TRUE; - for (i = 0; i < get_atomtype_ntypes(atype); i++) + for (int i = 0; i < gmx::ssize(atypes); i++) { - have_atomnumber = have_atomnumber && (get_atomtype_atomnumber(i, atype) >= 0); + have_atomnumber = have_atomnumber && (atypes.atomNumberFromAtomType(i) >= 0); } if (!have_atomnumber && ir->bQMMM) { @@ -1960,11 +1958,10 @@ int gmx_grompp(int argc, char *argv[]) setup_cmap(plist[F_CMAP].cmakeGridSpacing, plist[F_CMAP].cmapAngles, plist[F_CMAP].cmap, &sys.ffparams.cmap_grid); } - set_wall_atomtype(atype, opts, ir, wi); + set_wall_atomtype(&atypes, opts, ir, wi); if (bRenum) { - renum_atype(plist, &sys, ir->wall_atomtype, atype, bVerbose); - get_atomtype_ntypes(atype); + atypes.renumberTypes(plist, &sys, ir->wall_atomtype, bVerbose); } if (ir->implicit_solvent) @@ -1973,11 +1970,11 @@ int gmx_grompp(int argc, char *argv[]) } /* PELA: Copy the atomtype data to the topology atomtype list */ - copy_atomtype_atomtypes(atype, &(sys.atomtypes)); + atypes.copyTot_atomtypes(&(sys.atomtypes)); if (debug) { - pr_symtab(debug, 0, "After renum_atype", &sys.symtab); + pr_symtab(debug, 0, "After atype.renumberTypes", &sys.symtab); } if (bVerbose) @@ -1985,7 +1982,7 @@ int gmx_grompp(int argc, char *argv[]) fprintf(stderr, "converting bonded parameters...\n"); } - ntype = get_atomtype_ntypes(atype); + const int ntype = atypes.size(); convertInteractionTypeParameters(ntype, plist, mi, intermolecular_interactions.get(), comb, reppow, fudgeQQ, &sys); @@ -2329,7 +2326,6 @@ int gmx_grompp(int argc, char *argv[]) // fullCleanUp can't be called. mol.partialCleanUp(); } - done_atomtype(atype); done_inputrec_strings(); output_env_done(oenv); diff --git a/src/gromacs/gmxpreprocess/nm2type.cpp b/src/gromacs/gmxpreprocess/nm2type.cpp index 748e435b02..27c0a8a1c2 100644 --- a/src/gromacs/gmxpreprocess/nm2type.cpp +++ b/src/gromacs/gmxpreprocess/nm2type.cpp @@ -192,19 +192,18 @@ static int match_str(const char *atom, const char *template_string) } int nm2type(int nnm, t_nm2type nm2t[], t_symtab *tab, t_atoms *atoms, - gpp_atomtype *atype, int *nbonds, InteractionTypeParameters *bonds) + PreprocessingAtomTypes *atype, int *nbonds, InteractionTypeParameters *bonds) { int cur = 0; #define prev (1-cur) - int i, j, k, m, n, nresolved, nb, maxbond, ai, aj, best, im, nqual[2][ematchNR]; + int nresolved, nb, maxbond, nqual[2][ematchNR]; int *bbb, *n_mask, *m_mask, **match; - char *aname_i, *aname_m, *aname_n, *type; - double qq, mm; + char *aname_i, *aname_n; t_param *param; snew(param, 1); maxbond = 0; - for (i = 0; (i < atoms->nr); i++) + for (int i = 0; (i < atoms->nr); i++) { maxbond = std::max(maxbond, nbonds[i]); } @@ -216,20 +215,20 @@ int nm2type(int nnm, t_nm2type nm2t[], t_symtab *tab, t_atoms *atoms, snew(n_mask, maxbond); snew(m_mask, maxbond); snew(match, maxbond); - for (i = 0; (i < maxbond); i++) + for (int i = 0; (i < maxbond); i++) { snew(match[i], maxbond); } nresolved = 0; - for (i = 0; (i < atoms->nr); i++) + for (int i = 0; (i < atoms->nr); i++) { aname_i = *atoms->atomname[i]; nb = 0; - for (j = 0; (j < bonds->nr); j++) + for (int j = 0; (j < bonds->nr); j++) { - ai = bonds->param[j].ai(); - aj = bonds->param[j].aj(); + int ai = bonds->param[j].ai(); + int aj = bonds->param[j].aj(); if (ai == i) { bbb[nb++] = aj; @@ -247,52 +246,52 @@ int nm2type(int nnm, t_nm2type nm2t[], t_symtab *tab, t_atoms *atoms, if (debug) { fprintf(debug, "%4s has bonds to", aname_i); - for (j = 0; (j < nb); j++) + for (int j = 0; (j < nb); j++) { fprintf(debug, " %4s", *atoms->atomname[bbb[j]]); } fprintf(debug, "\n"); } - best = -1; - for (k = 0; (k < ematchNR); k++) + int best = -1; + for (int k = 0; (k < ematchNR); k++) { nqual[prev][k] = 0; } /* First check for names */ - for (k = 0; (k < nnm); k++) + for (int k = 0; (k < nnm); k++) { if (nm2t[k].nbonds == nb) { - im = match_str(*atoms->atomname[i], nm2t[k].elem); + int im = match_str(*atoms->atomname[i], nm2t[k].elem); if (im > ematchWild) { - for (j = 0; (j < ematchNR); j++) + for (int j = 0; (j < ematchNR); j++) { nqual[cur][j] = 0; } /* Fill a matrix with matching quality */ - for (m = 0; (m < nb); m++) + for (int m = 0; (m < nb); m++) { - aname_m = *atoms->atomname[bbb[m]]; - for (n = 0; (n < nb); n++) + const char *aname_m = *atoms->atomname[bbb[m]]; + for (int n = 0; (n < nb); n++) { aname_n = nm2t[k].bond[n]; match[m][n] = match_str(aname_m, aname_n); } } /* Now pick the best matches */ - for (m = 0; (m < nb); m++) + for (int m = 0; (m < nb); m++) { n_mask[m] = 0; m_mask[m] = 0; } - for (j = ematchNR-1; (j > 0); j--) + for (int j = ematchNR-1; (j > 0); j--) { - for (m = 0; (m < nb); m++) + for (int m = 0; (m < nb); m++) { - for (n = 0; (n < nb); n++) + for (int n = 0; (n < nb); n++) { if ((n_mask[n] == 0) && (m_mask[m] == 0) && @@ -327,19 +326,20 @@ int nm2type(int nnm, t_nm2type nm2t[], t_symtab *tab, t_atoms *atoms, } if (best != -1) { - int atomnr = 0; - real alpha = 0; + int atomnr = 0; + real alpha = 0; - qq = nm2t[best].q; - mm = nm2t[best].m; - type = nm2t[best].type; + double qq = nm2t[best].q; + double mm = nm2t[best].m; + const char *type = nm2t[best].type; - if ((k = get_atomtype_type(type, atype)) == NOTSET) + int k; + if ((k = atype->atomTypeFromName(type)) == NOTSET) { atoms->atom[i].qB = alpha; atoms->atom[i].m = atoms->atom[i].mB = mm; - k = add_atomtype(atype, tab, &(atoms->atom[i]), type, param, - atoms->atom[i].type, atomnr); + k = atype->addType(tab, atoms->atom[i], type, param, + atoms->atom[i].type, atomnr); } atoms->atom[i].type = k; atoms->atom[i].typeB = k; diff --git a/src/gromacs/gmxpreprocess/nm2type.h b/src/gromacs/gmxpreprocess/nm2type.h index 5d6fa3d0b4..f982ececb7 100644 --- a/src/gromacs/gmxpreprocess/nm2type.h +++ b/src/gromacs/gmxpreprocess/nm2type.h @@ -39,7 +39,7 @@ #include -struct gpp_atomtype; +class PreprocessingAtomTypes; struct t_atoms; struct InteractionTypeParameters; struct t_symtab; @@ -62,7 +62,7 @@ void dump_nm2type(FILE *fp, int nnm, t_nm2type nm2t[]); /* Dump the database for debugging. Can be reread by the program */ int nm2type(int nnm, t_nm2type nm2t[], t_symtab *tab, t_atoms *atoms, - gpp_atomtype *atype, int *nbonds, InteractionTypeParameters *bond); + PreprocessingAtomTypes *atype, int *nbonds, InteractionTypeParameters *bond); /* Try to determine the atomtype (force field dependent) for the atoms * with help of the bond list */ diff --git a/src/gromacs/gmxpreprocess/pdb2gmx.cpp b/src/gromacs/gmxpreprocess/pdb2gmx.cpp index aaea0e566d..ca33d56ca8 100644 --- a/src/gromacs/gmxpreprocess/pdb2gmx.cpp +++ b/src/gromacs/gmxpreprocess/pdb2gmx.cpp @@ -1912,7 +1912,7 @@ int pdb2gmx::run() check_occupancy(&pdba_all, inputConfFile_.c_str(), bVerbose_); /* Read atomtypes... */ - gpp_atomtype *atype = read_atype(ffdir_, &symtab); + PreprocessingAtomTypes atype = read_atype(ffdir_, &symtab); /* read residue database */ printf("Reading residue database... (%s)\n", forcefield_); @@ -1920,7 +1920,7 @@ int pdb2gmx::run() std::vector rtpFFDB; for (const auto &filename : rtpf) { - readResidueDatabase(filename, &rtpFFDB, atype, &symtab, false); + readResidueDatabase(filename, &rtpFFDB, &atype, &symtab, false); } if (bNewRTP_) { @@ -1938,8 +1938,8 @@ int pdb2gmx::run() std::vector ntdb; std::vector ctdb; std::vector tdblist; - int nNtdb = read_ter_db(ffdir_, 'n', &ntdb, atype); - int nCtdb = read_ter_db(ffdir_, 'c', &ctdb, atype); + int nNtdb = read_ter_db(ffdir_, 'n', &ntdb, &atype); + int nCtdb = read_ter_db(ffdir_, 'c', &ctdb, &atype); FILE *top_file = gmx_fio_fopen(topologyFile_.c_str(), "w"); @@ -2249,7 +2249,7 @@ int pdb2gmx::run() top_file2 = top_file; } - pdb2top(top_file2, posre_fn.c_str(), molname.c_str(), pdba, &x, atype, &symtab, + pdb2top(top_file2, posre_fn.c_str(), molname.c_str(), pdba, &x, &atype, &symtab, rtpFFDB, restp_chain, hb_chain, bAllowMissing_, @@ -2272,7 +2272,6 @@ int pdb2gmx::run() cc->pdba = pdba; cc->x = x; } - done_atomtype(atype); if (watermodel_ == nullptr) { diff --git a/src/gromacs/gmxpreprocess/pdb2top.cpp b/src/gromacs/gmxpreprocess/pdb2top.cpp index ed6a9c7d4d..f2f24f03d4 100644 --- a/src/gromacs/gmxpreprocess/pdb2top.cpp +++ b/src/gromacs/gmxpreprocess/pdb2top.cpp @@ -673,7 +673,7 @@ void print_top_mols(FILE *out, void write_top(FILE *out, const char *pr, const char *molname, t_atoms *at, bool bRTPresname, int bts[], gmx::ArrayRef plist, t_excls excls[], - gpp_atomtype *atype, int *cgnr, int nrexcl) + PreprocessingAtomTypes *atype, int *cgnr, int nrexcl) /* NOTE: nrexcl is not the size of *excl! */ { if (at && atype && cgnr) @@ -1465,7 +1465,7 @@ scrub_charge_groups(int *cgnr, int natoms) void pdb2top(FILE *top_file, const char *posre_fn, const char *molname, t_atoms *atoms, - std::vector *x, gpp_atomtype *atype, t_symtab *tab, + std::vector *x, PreprocessingAtomTypes *atype, t_symtab *tab, gmx::ArrayRef rtpFFDB, gmx::ArrayRef usedPpResidues, gmx::ArrayRef globalPatches, diff --git a/src/gromacs/gmxpreprocess/pdb2top.h b/src/gromacs/gmxpreprocess/pdb2top.h index b9c625c3a4..072371e15e 100644 --- a/src/gromacs/gmxpreprocess/pdb2top.h +++ b/src/gromacs/gmxpreprocess/pdb2top.h @@ -45,7 +45,7 @@ #include "gromacs/math/vectypes.h" #include "gromacs/utility/arrayref.h" -struct gpp_atomtype; +class PreprocessingAtomTypes; struct t_atoms; struct t_excls; struct MoleculePatchDatabase; @@ -114,13 +114,13 @@ void print_top_mols(FILE *out, void write_top(FILE *out, const char *pr, const char *molname, t_atoms *at, bool bRTPresname, int bts[], gmx::ArrayRef plist, t_excls excls[], - gpp_atomtype *atype, int *cgnr, int nrexcl); + PreprocessingAtomTypes *atype, int *cgnr, int nrexcl); /* NOTE: nrexcl is not the size of *excl! */ void pdb2top(FILE *top_file, const char *posre_fn, const char *molname, t_atoms *atoms, std::vector *x, - gpp_atomtype *atype, t_symtab *tab, + PreprocessingAtomTypes *atype, t_symtab *tab, gmx::ArrayRef rtpFFDB, gmx::ArrayRef usedPpResidues, gmx::ArrayRef globalPatches, diff --git a/src/gromacs/gmxpreprocess/resall.cpp b/src/gromacs/gmxpreprocess/resall.cpp index 1f4ad212de..37e6024d6f 100644 --- a/src/gromacs/gmxpreprocess/resall.cpp +++ b/src/gromacs/gmxpreprocess/resall.cpp @@ -61,20 +61,19 @@ #include "hackblock.h" -gpp_atomtype *read_atype(const char *ffdir, t_symtab *tab) +PreprocessingAtomTypes read_atype(const char *ffdir, t_symtab *tab) { FILE *in; char buf[STRLEN], name[STRLEN]; double m; int nratt = 0; - gpp_atomtype *at; t_atom *a; t_param *nb; std::vector files = fflib_search_file_end(ffdir, ".atp", TRUE); - at = init_atomtype(); snew(a, 1); snew(nb, 1); + PreprocessingAtomTypes at; for (const auto &filename : files) { @@ -95,7 +94,7 @@ gpp_atomtype *read_atype(const char *ffdir, t_symtab *tab) if (sscanf(buf, "%s%lf", name, &m) == 2) { a->m = m; - add_atomtype(at, tab, a, name, nb, 0, 0); + at.addType(tab, *a, name, nb, 0, 0); fprintf(stderr, "\rAtomtype %d", ++nratt); fflush(stderr); } @@ -108,12 +107,10 @@ gpp_atomtype *read_atype(const char *ffdir, t_symtab *tab) } fprintf(stderr, "\n"); sfree(a); - sfree(nb); - return at; } -static void print_resatoms(FILE *out, gpp_atomtype *atype, const PreprocessResidue &rtpDBEntry) +static void print_resatoms(FILE *out, const PreprocessingAtomTypes &atype, const PreprocessResidue &rtpDBEntry) { fprintf(out, "[ %s ]\n", rtpDBEntry.resname.c_str()); fprintf(out, " [ atoms ]\n"); @@ -121,7 +118,7 @@ static void print_resatoms(FILE *out, gpp_atomtype *atype, const PreprocessResid for (int j = 0; (j < rtpDBEntry.natom()); j++) { int tp = rtpDBEntry.atom[j].type; - const char *tpnm = get_atomtype_name(tp, atype); + const char *tpnm = atype.atomNameFromAtomType(tp); if (tpnm == nullptr) { gmx_fatal(FARGS, "Incorrect atomtype (%d)", tp); @@ -132,7 +129,7 @@ static void print_resatoms(FILE *out, gpp_atomtype *atype, const PreprocessResid } static bool read_atoms(FILE *in, char *line, - PreprocessResidue *r0, t_symtab *tab, gpp_atomtype *atype) + PreprocessResidue *r0, t_symtab *tab, PreprocessingAtomTypes *atype) { int cg; char buf[256], buf1[256]; @@ -152,14 +149,14 @@ static bool read_atoms(FILE *in, char *line, r0->atom.emplace_back(); r0->atom.back().q = q; r0->cgnr.push_back(cg); - int j = get_atomtype_type(buf1, atype); + int j = atype->atomTypeFromName(buf1); if (j == NOTSET) { gmx_fatal(FARGS, "Atom type %s (residue %s) not found in atomtype " "database", buf1, r0->resname.c_str()); } r0->atom.back().type = j; - r0->atom.back().m = get_atomtype_massA(j, atype); + r0->atom.back().m = atype->atomMassAFromAtomType(j); } return TRUE; @@ -264,7 +261,7 @@ static void print_resall_header(FILE *out, gmx::ArrayRef rtpDBEntry, - gpp_atomtype *atype) + const PreprocessingAtomTypes &atype) { if (rtpDBEntry.empty()) { @@ -287,7 +284,7 @@ void print_resall(FILE *out, gmx::ArrayRef rtpDBEntry, } void readResidueDatabase(const std::string &rrdb, std::vector *rtpDBEntry, - gpp_atomtype *atype, t_symtab *tab, + PreprocessingAtomTypes *atype, t_symtab *tab, bool bAllowOverrideRTP) { FILE *in; diff --git a/src/gromacs/gmxpreprocess/resall.h b/src/gromacs/gmxpreprocess/resall.h index 7b623b9325..c5022ca246 100644 --- a/src/gromacs/gmxpreprocess/resall.h +++ b/src/gromacs/gmxpreprocess/resall.h @@ -47,7 +47,7 @@ #include "gromacs/utility/arrayref.h" -struct gpp_atomtype; +class PreprocessingAtomTypes; struct PreprocessResidue; struct t_symtab; @@ -81,22 +81,22 @@ getDatabaseEntry(const std::string &rtpname, gmx::ArrayRef *rtpDBEntry, - gpp_atomtype *atype, - t_symtab *tab, - bool bAllowOverrideRTP); +void readResidueDatabase(const std::string &resdb, + std::vector *rtpDBEntry, + PreprocessingAtomTypes *atype, + t_symtab *tab, + bool bAllowOverrideRTP); /*! \brief * Print out database. @@ -106,6 +106,6 @@ void readResidueDatabase(const std::string &resdb, * \param[in] atype Atom type information. */ void print_resall(FILE *out, gmx::ArrayRef rtpDBEntry, - gpp_atomtype *atype); + const PreprocessingAtomTypes &atype); #endif diff --git a/src/gromacs/gmxpreprocess/ter_db.cpp b/src/gromacs/gmxpreprocess/ter_db.cpp index 8989d23543..9677b742ce 100644 --- a/src/gromacs/gmxpreprocess/ter_db.cpp +++ b/src/gromacs/gmxpreprocess/ter_db.cpp @@ -95,7 +95,7 @@ static int find_kw(char *keyw) #define FATAL() gmx_fatal(FARGS, "Reading Termini Database: not enough items on line\n%s", line) static void read_atom(char *line, bool bAdd, - std::string *nname, t_atom *a, gpp_atomtype *atype, int *cgnr) + std::string *nname, t_atom *a, PreprocessingAtomTypes *atype, int *cgnr) { int nr, i; char buf[5][30]; @@ -133,7 +133,7 @@ static void read_atom(char *line, bool bAdd, *nname = ""; } } - a->type = get_atomtype_type(buf[i++], atype); + a->type = atype->atomTypeFromName(buf[i++]); sscanf(buf[i++], "%lf", &m); a->m = m; sscanf(buf[i++], "%lf", &q); @@ -148,14 +148,14 @@ static void read_atom(char *line, bool bAdd, } } -static void print_atom(FILE *out, const t_atom *a, gpp_atomtype *atype) +static void print_atom(FILE *out, const t_atom &a, PreprocessingAtomTypes *atype) { fprintf(out, "\t%s\t%g\t%g\n", - get_atomtype_name(a->type, atype), a->m, a->q); + atype->atomNameFromAtomType(a.type), a.m, a.q); } static void print_ter_db(const char *ff, char C, gmx::ArrayRef tb, - gpp_atomtype *atype) + PreprocessingAtomTypes *atype) { std::string buf = gmx::formatString("%s-%c.tdb", ff, C); FILE *out = gmx_fio_fopen(buf.c_str(), "w"); @@ -174,7 +174,7 @@ static void print_ter_db(const char *ff, char C, gmx::ArrayRef *tbptr, - gpp_atomtype *atype) +static void read_ter_db_file(const char *fn, + std::vector *tbptr, + PreprocessingAtomTypes *atype) { char filebase[STRLEN]; char header[STRLEN], buf[STRLEN], line[STRLEN]; @@ -359,7 +359,7 @@ static void read_ter_db_file(const char * } int read_ter_db(const char *ffdir, char ter, - std::vector *tbptr, gpp_atomtype *atype) + std::vector *tbptr, PreprocessingAtomTypes *atype) { std::string ext = gmx::formatString(".%c.tdb", ter); diff --git a/src/gromacs/gmxpreprocess/ter_db.h b/src/gromacs/gmxpreprocess/ter_db.h index 9961c77641..c3e531e7fd 100644 --- a/src/gromacs/gmxpreprocess/ter_db.h +++ b/src/gromacs/gmxpreprocess/ter_db.h @@ -42,7 +42,7 @@ #include "gromacs/utility/arrayref.h" -struct gpp_atomtype; +class PreprocessingAtomTypes; struct MoleculePatchDatabase; /*! \brief @@ -55,7 +55,7 @@ struct MoleculePatchDatabase; * \returns Number of entries entered into database. */ int read_ter_db(const char *ffdir, char ter, - std::vector *tbptr, gpp_atomtype *atype); + std::vector *tbptr, PreprocessingAtomTypes *atype); /*! \brief * Return entries for modification blocks that match a residue name. diff --git a/src/gromacs/gmxpreprocess/tomorse.cpp b/src/gromacs/gmxpreprocess/tomorse.cpp index 9bd2c09262..934a773b0c 100644 --- a/src/gromacs/gmxpreprocess/tomorse.cpp +++ b/src/gromacs/gmxpreprocess/tomorse.cpp @@ -99,7 +99,7 @@ static t_2morse *read_dissociation_energies(int *n2morse) return t2m; } -static int nequal(char *a1, char *a2) +static int nequal(const char *a1, const char *a2) { int i; @@ -122,7 +122,7 @@ static int nequal(char *a1, char *a2) return i; } -static real search_e_diss(int n2m, t_2morse t2m[], char *ai, char *aj) +static real search_e_diss(int n2m, t_2morse t2m[], const char *ai, const char *aj) { int i; int ibest = -1; @@ -187,7 +187,7 @@ static real search_e_diss(int n2m, t_2morse t2m[], char *ai, char *aj) } } -void convert_harmonics(gmx::ArrayRef mols, gpp_atomtype *atype) +void convert_harmonics(gmx::ArrayRef mols, PreprocessingAtomTypes *atype) { int n2m; t_2morse *t2m; @@ -226,8 +226,8 @@ void convert_harmonics(gmx::ArrayRef mols, gpp_atomtype *at int nj = mol.plist[bb].param[j].aj(); real edis = search_e_diss(n2m, t2m, - get_atomtype_name(mol.atoms.atom[ni].type, atype), - get_atomtype_name(mol.atoms.atom[nj].type, atype)); + atype->atomNameFromAtomType(mol.atoms.atom[ni].type), + atype->atomNameFromAtomType(mol.atoms.atom[nj].type)); if (edis != 0) { bRemoveHarm[j] = true; diff --git a/src/gromacs/gmxpreprocess/tomorse.h b/src/gromacs/gmxpreprocess/tomorse.h index 9d2f510493..eb20491064 100644 --- a/src/gromacs/gmxpreprocess/tomorse.h +++ b/src/gromacs/gmxpreprocess/tomorse.h @@ -40,9 +40,9 @@ #include "gromacs/utility/arrayref.h" -struct gpp_atomtype; +class PreprocessingAtomTypes; struct MoleculeInformation; -void convert_harmonics(gmx::ArrayRef mols, gpp_atomtype *atype); +void convert_harmonics(gmx::ArrayRef mols, PreprocessingAtomTypes *atype); #endif diff --git a/src/gromacs/gmxpreprocess/topio.cpp b/src/gromacs/gmxpreprocess/topio.cpp index ae0f82e37d..eb0beb201f 100644 --- a/src/gromacs/gmxpreprocess/topio.cpp +++ b/src/gromacs/gmxpreprocess/topio.cpp @@ -85,16 +85,15 @@ #define OPENDIR '[' /* starting sign for directive */ #define CLOSEDIR ']' /* ending sign for directive */ -static void gen_pairs(InteractionTypeParameters *nbs, InteractionTypeParameters *pairs, real fudge, int comb) +static void gen_pairs(const InteractionTypeParameters &nbs, InteractionTypeParameters *pairs, real fudge, int comb) { - int i, j, ntp, nrfp, nrfpA, nrfpB, nnn; real scaling; - ntp = nbs->nr; - nnn = static_cast(std::sqrt(static_cast(ntp))); + int ntp = nbs.nr; + int nnn = static_cast(std::sqrt(static_cast(ntp))); GMX_ASSERT(nnn * nnn == ntp, "Number of pairs of generated non-bonded parameters should be a perfect square"); - nrfp = NRFP(F_LJ); - nrfpA = interaction_function[F_LJ14].nrfpA; - nrfpB = interaction_function[F_LJ14].nrfpB; + int nrfp = NRFP(F_LJ); + int nrfpA = interaction_function[F_LJ14].nrfpA; + int nrfpB = interaction_function[F_LJ14].nrfpB; pairs->nr = ntp; if ((nrfp != nrfpA) || (nrfpA != nrfpB)) @@ -104,7 +103,7 @@ static void gen_pairs(InteractionTypeParameters *nbs, InteractionTypeParameters fprintf(stderr, "Generating 1-4 interactions: fudge = %g\n", fudge); snew(pairs->param, pairs->nr); - for (i = 0; (i < ntp); i++) + for (int i = 0; (i < ntp); i++) { /* Copy param.a */ pairs->param[i].a[0] = i / nnn; @@ -113,7 +112,7 @@ static void gen_pairs(InteractionTypeParameters *nbs, InteractionTypeParameters - for (j = 0; (j < nrfp); j++) + for (int j = 0; (j < nrfp); j++) { /* If we are using sigma/epsilon values, only the epsilon values * should be scaled, but not sigma. @@ -128,12 +127,12 @@ static void gen_pairs(InteractionTypeParameters *nbs, InteractionTypeParameters scaling = fudge; } - pairs->param[i].c[j] = scaling*nbs->param[i].c[j]; + pairs->param[i].c[j] = scaling*nbs.param[i].c[j]; /* NOTE: this should be clear to the compiler, but some gcc 5.2 versions * issue false positive warnings for the pairs->param.c[] indexing below. */ assert(2*nrfp <= MAXFORCEPARAM); - pairs->param[i].c[nrfp+j] = scaling*nbs->param[i].c[j]; + pairs->param[i].c[nrfp+j] = scaling*nbs.param[i].c[j]; } } } @@ -377,8 +376,7 @@ static void make_atoms_sys(gmx::ArrayRef molblock, static char **read_topol(const char *infile, const char *outfile, const char *define, const char *include, t_symtab *symtab, - gpp_atomtype *atype, - int *nrmols, + PreprocessingAtomTypes *atypes, std::vector *molinfo, std::unique_ptr *intermolecular_interactions, gmx::ArrayRef plist, @@ -399,7 +397,7 @@ static char **read_topol(const char *infile, const char *outfile, char line[STRLEN], errbuf[256], comb_str[256], nb_str[256]; char genpairs[32]; char *dirstr, *dummy2; - int nrcopies, nmol, nscan, ncombs, ncopy; + int nrcopies, nscan, ncombs, ncopy; double fLJ, fQQ, fPOW; MoleculeInformation *mi0 = nullptr; DirStack *DS; @@ -443,7 +441,6 @@ static char **read_topol(const char *infile, const char *outfile, /* some local variables */ DS_Init(&DS); /* directive stack */ - nmol = 0; /* no molecules yet... */ d = Directive::d_invalid; /* first thing should be a directive */ nbparam = nullptr; /* The temporary non-bonded matrix */ pair = nullptr; /* The temporary pair interaction matrix */ @@ -639,7 +636,7 @@ static char **read_topol(const char *infile, const char *outfile, break; case Directive::d_atomtypes: - push_at(symtab, atype, batype, pline, nb_funct, + push_at(symtab, atypes, batype, pline, nb_funct, &nbparam, bGenPairs ? &pair : nullptr, wi); break; @@ -652,11 +649,11 @@ static char **read_topol(const char *infile, const char *outfile, case Directive::d_pairtypes: if (bGenPairs) { - push_nbt(d, pair, atype, pline, F_LJ14, wi); + push_nbt(d, pair, atypes, pline, F_LJ14, wi); } else { - push_bt(d, plist, 2, atype, nullptr, pline, wi); + push_bt(d, plist, 2, atypes, nullptr, pline, wi); } break; case Directive::d_angletypes: @@ -668,7 +665,7 @@ static char **read_topol(const char *infile, const char *outfile, break; case Directive::d_nonbond_params: - push_nbt(d, nbparam, atype, pline, nb_funct, wi); + push_nbt(d, nbparam, atypes, pline, nb_funct, wi); break; case Directive::d_implicit_genborn_params: @@ -684,7 +681,7 @@ static char **read_topol(const char *infile, const char *outfile, break; case Directive::d_cmaptypes: - push_cmaptype(d, plist, 5, atype, batype, pline, wi); + push_cmaptype(d, plist, 5, atypes, batype, pline, wi); break; case Directive::d_moleculetype: @@ -698,19 +695,19 @@ static char **read_topol(const char *infile, const char *outfile, opts->couple_lam1 == ecouplamNONE || opts->couple_lam1 == ecouplamQ)) { - dcatt = add_atomtype_decoupled(symtab, atype, + dcatt = add_atomtype_decoupled(symtab, atypes, &nbparam, bGenPairs ? &pair : nullptr); } - ntype = get_atomtype_ntypes(atype); + ntype = atypes->size(); ncombs = (ntype*(ntype+1))/2; - generate_nbparams(*combination_rule, nb_funct, &(plist[nb_funct]), atype, wi); + generate_nbparams(*combination_rule, nb_funct, &(plist[nb_funct]), atypes, wi); ncopy = copy_nbparams(nbparam, nb_funct, &(plist[nb_funct]), ntype); fprintf(stderr, "Generated %d of the %d non-bonded parameter combinations\n", ncombs-ncopy, ncombs); free_nbparam(nbparam, ntype); if (bGenPairs) { - gen_pairs(&(plist[nb_funct]), &(plist[F_LJ14]), fudgeLJ, *combination_rule); + gen_pairs((plist[nb_funct]), &(plist[F_LJ14]), fudgeLJ, *combination_rule); ncopy = copy_nbparams(pair, nb_funct, &(plist[F_LJ14]), ntype); fprintf(stderr, "Generated %d of the %d 1-4 parameter combinations\n", ncombs-ncopy, ncombs); @@ -722,7 +719,6 @@ static char **read_topol(const char *infile, const char *outfile, } push_molt(symtab, molinfo, pline, wi); - nmol = molinfo->size(); exclusionBlocks.emplace_back(); mi0 = &molinfo->back(); mi0->atoms.haveMass = TRUE; @@ -733,17 +729,17 @@ static char **read_topol(const char *infile, const char *outfile, break; } case Directive::d_atoms: - push_atom(symtab, &(mi0->cgs), &(mi0->atoms), atype, pline, &lastcg, wi); + push_atom(symtab, &(mi0->cgs), &(mi0->atoms), atypes, pline, &lastcg, wi); break; case Directive::d_pairs: GMX_RELEASE_ASSERT(mi0, "Need to have a valid MoleculeInformation object to work on"); - push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, FALSE, + push_bond(d, plist, mi0->plist, &(mi0->atoms), atypes, pline, FALSE, bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi); break; case Directive::d_pairs_nb: GMX_RELEASE_ASSERT(mi0, "Need to have a valid MoleculeInformation object to work on"); - push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, FALSE, + push_bond(d, plist, mi0->plist, &(mi0->atoms), atypes, pline, FALSE, FALSE, 1.0, bZero, &bWarn_copy_A_B, wi); break; @@ -765,12 +761,12 @@ static char **read_topol(const char *infile, const char *outfile, case Directive::d_water_polarization: case Directive::d_thole_polarization: GMX_RELEASE_ASSERT(mi0, "Need to have a valid MoleculeInformation object to work on"); - push_bond(d, plist, mi0->plist, &(mi0->atoms), atype, pline, TRUE, + push_bond(d, plist, mi0->plist, &(mi0->atoms), atypes, pline, TRUE, bGenPairs, *fudgeQQ, bZero, &bWarn_copy_A_B, wi); break; case Directive::d_cmap: GMX_RELEASE_ASSERT(mi0, "Need to have a valid MoleculeInformation object to work on"); - push_cmap(d, plist, mi0->plist, &(mi0->atoms), atype, pline, wi); + push_cmap(d, plist, mi0->plist, &(mi0->atoms), atypes, pline, wi); break; case Directive::d_vsitesn: @@ -933,29 +929,26 @@ static char **read_topol(const char *infile, const char *outfile, sfree(intermolecular_interactions->get()->atoms.atom); } - *nrmols = nmol; - return title; } -char **do_top(bool bVerbose, - const char *topfile, - const char *topppfile, - t_gromppopts *opts, - bool bZero, - t_symtab *symtab, - gmx::ArrayRef plist, - int *combination_rule, - double *repulsion_power, - real *fudgeQQ, - gpp_atomtype *atype, - int *nrmols, - std::vector *molinfo, - std::unique_ptr *intermolecular_interactions, - const t_inputrec *ir, - std::vector *molblock, - bool *ffParametrizedWithHBondConstraints, - warninp *wi) +char **do_top(bool bVerbose, + const char *topfile, + const char *topppfile, + t_gromppopts *opts, + bool bZero, + t_symtab *symtab, + gmx::ArrayRef plist, + int *combination_rule, + double *repulsion_power, + real *fudgeQQ, + PreprocessingAtomTypes *atypes, + std::vector *molinfo, + std::unique_ptr *intermolecular_interactions, + const t_inputrec *ir, + std::vector *molblock, + bool *ffParametrizedWithHBondConstraints, + warninp *wi) { /* Tmpfile might contain a long path */ const char *tmpfile; @@ -975,8 +968,8 @@ char **do_top(bool bVerbose, printf("processing topology...\n"); } title = read_topol(topfile, tmpfile, opts->define, opts->include, - symtab, atype, - nrmols, molinfo, intermolecular_interactions, + symtab, atypes, + molinfo, intermolecular_interactions, plist, combination_rule, repulsion_power, opts, fudgeQQ, molblock, ffParametrizedWithHBondConstraints, diff --git a/src/gromacs/gmxpreprocess/topio.h b/src/gromacs/gmxpreprocess/topio.h index 746aafe3a2..c5e5386e3a 100644 --- a/src/gromacs/gmxpreprocess/topio.h +++ b/src/gromacs/gmxpreprocess/topio.h @@ -46,7 +46,7 @@ struct gmx_molblock_t; struct gmx_mtop_t; -struct gpp_atomtype; +class PreprocessingAtomTypes; struct t_gromppopts; struct t_inputrec; struct MoleculeInformation; @@ -59,24 +59,23 @@ 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, - gmx::ArrayRef plist, - int *combination_rule, - double *repulsion_power, - real *fudgeQQ, - gpp_atomtype *atype, - int *nrmols, - std::vector *molinfo, - std::unique_ptr *intermolecular_interactions, - const t_inputrec *ir, - std::vector *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, + gmx::ArrayRef plist, + int *combination_rule, + double *repulsion_power, + real *fudgeQQ, + PreprocessingAtomTypes *atype, + std::vector *molinfo, + std::unique_ptr *intermolecular_interactions, + const t_inputrec *ir, + std::vector *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); diff --git a/src/gromacs/gmxpreprocess/toppush.cpp b/src/gromacs/gmxpreprocess/toppush.cpp index ca06defb24..0d6dbcd7cc 100644 --- a/src/gromacs/gmxpreprocess/toppush.cpp +++ b/src/gromacs/gmxpreprocess/toppush.cpp @@ -66,15 +66,18 @@ #include "gromacs/utility/smalloc.h" #include "gromacs/utility/stringutil.h" -void generate_nbparams(int comb, int ftype, InteractionTypeParameters *plist, gpp_atomtype *atype, - warninp *wi) +void generate_nbparams(int comb, + int ftype, + InteractionTypeParameters *plist, + PreprocessingAtomTypes *atypes, + warninp *wi) { - int i, j, k = -1, nf; + int k = -1; int nr, nrfp; real c, bi, bj, ci, cj, ci0, ci1, ci2, cj0, cj1, cj2; /* Lean mean shortcuts */ - nr = get_atomtype_ntypes(atype); + nr = atypes->size(); nrfp = NRFP(ftype); snew(plist->param, nr*nr); plist->nr = nr*nr; @@ -87,14 +90,14 @@ void generate_nbparams(int comb, int ftype, InteractionTypeParameters *plist, gp { case eCOMB_GEOMETRIC: /* Gromos rules */ - for (i = k = 0; (i < nr); i++) + for (int i = k = 0; (i < nr); i++) { - for (j = 0; (j < nr); j++, k++) + for (int j = 0; (j < nr); j++, k++) { - for (nf = 0; (nf < nrfp); nf++) + for (int nf = 0; (nf < nrfp); nf++) { - ci = get_atomtype_nbparam(i, nf, atype); - cj = get_atomtype_nbparam(j, nf, atype); + ci = atypes->atomNonBondedParamFromAtomType(i, nf); + cj = atypes->atomNonBondedParamFromAtomType(j, nf); c = std::sqrt(ci * cj); plist->param[k].c[nf] = c; } @@ -104,14 +107,14 @@ void generate_nbparams(int comb, int ftype, InteractionTypeParameters *plist, gp case eCOMB_ARITHMETIC: /* c0 and c1 are sigma and epsilon */ - for (i = k = 0; (i < nr); i++) + for (int i = k = 0; (i < nr); i++) { - for (j = 0; (j < nr); j++, k++) + for (int j = 0; (j < nr); j++, k++) { - ci0 = get_atomtype_nbparam(i, 0, atype); - cj0 = get_atomtype_nbparam(j, 0, atype); - ci1 = get_atomtype_nbparam(i, 1, atype); - cj1 = get_atomtype_nbparam(j, 1, atype); + ci0 = atypes->atomNonBondedParamFromAtomType(i, 0); + cj0 = atypes->atomNonBondedParamFromAtomType(j, 0); + ci1 = atypes->atomNonBondedParamFromAtomType(i, 1); + cj1 = atypes->atomNonBondedParamFromAtomType(j, 1); plist->param[k].c[0] = (fabs(ci0) + fabs(cj0))*0.5; /* Negative sigma signals that c6 should be set to zero later, * so we need to propagate that through the combination rules. @@ -127,14 +130,14 @@ void generate_nbparams(int comb, int ftype, InteractionTypeParameters *plist, gp break; case eCOMB_GEOM_SIG_EPS: /* c0 and c1 are sigma and epsilon */ - for (i = k = 0; (i < nr); i++) + for (int i = k = 0; (i < nr); i++) { - for (j = 0; (j < nr); j++, k++) + for (int j = 0; (j < nr); j++, k++) { - ci0 = get_atomtype_nbparam(i, 0, atype); - cj0 = get_atomtype_nbparam(j, 0, atype); - ci1 = get_atomtype_nbparam(i, 1, atype); - cj1 = get_atomtype_nbparam(j, 1, atype); + ci0 = atypes->atomNonBondedParamFromAtomType(i, 0); + cj0 = atypes->atomNonBondedParamFromAtomType(j, 0); + ci1 = atypes->atomNonBondedParamFromAtomType(i, 1); + cj1 = atypes->atomNonBondedParamFromAtomType(j, 1); plist->param[k].c[0] = std::sqrt(std::fabs(ci0*cj0)); /* Negative sigma signals that c6 should be set to zero later, * so we need to propagate that through the combination rules. @@ -160,16 +163,16 @@ void generate_nbparams(int comb, int ftype, InteractionTypeParameters *plist, gp case F_BHAM: /* Buckingham rules */ - for (i = k = 0; (i < nr); i++) + for (int i = k = 0; (i < nr); i++) { - for (j = 0; (j < nr); j++, k++) + for (int j = 0; (j < nr); j++, k++) { - ci0 = get_atomtype_nbparam(i, 0, atype); - cj0 = get_atomtype_nbparam(j, 0, atype); - ci2 = get_atomtype_nbparam(i, 2, atype); - cj2 = get_atomtype_nbparam(j, 2, atype); - bi = get_atomtype_nbparam(i, 1, atype); - bj = get_atomtype_nbparam(j, 1, atype); + ci0 = atypes->atomNonBondedParamFromAtomType(i, 0); + cj0 = atypes->atomNonBondedParamFromAtomType(j, 0); + ci2 = atypes->atomNonBondedParamFromAtomType(i, 2); + cj2 = atypes->atomNonBondedParamFromAtomType(j, 2); + bi = atypes->atomNonBondedParamFromAtomType(i, 1); + bj = atypes->atomNonBondedParamFromAtomType(j, 1); plist->param[k].c[0] = std::sqrt(ci0 * cj0); if ((bi == 0) || (bj == 0)) { @@ -201,11 +204,11 @@ struct t_nbparam real c[4]; }; -static void realloc_nb_params(gpp_atomtype *at, +static void realloc_nb_params(PreprocessingAtomTypes *atypes, t_nbparam ***nbparam, t_nbparam ***pair) { /* Add space in the non-bonded parameters matrix */ - int atnr = get_atomtype_ntypes(at); + int atnr = atypes->size(); srenew(*nbparam, atnr); snew((*nbparam)[atnr-1], atnr); if (pair) @@ -270,7 +273,7 @@ static void copy_B_from_A(int ftype, double *c) } } -void push_at (t_symtab *symtab, gpp_atomtype *at, gpp_bond_atomtype *bat, +void push_at (t_symtab *symtab, PreprocessingAtomTypes *at, gpp_bond_atomtype *bat, char *line, int nb_funct, t_nbparam ***nbparam, t_nbparam ***pair, warninp *wi) @@ -499,7 +502,7 @@ void push_at (t_symtab *symtab, gpp_atomtype *at, gpp_bond_atomtype *bat, auto message = gmx::formatString("Invalid function type %d in push_at", nb_funct); warning_error_and_exit(wi, message, FARGS); } - for (j = nfp0; (j < MAXFORCEPARAM); j++) + for (int j = nfp0; (j < MAXFORCEPARAM); j++) { c[j] = 0.0; } @@ -547,19 +550,19 @@ void push_at (t_symtab *symtab, gpp_atomtype *at, gpp_bond_atomtype *bat, } batype_nr = get_bond_atomtype_type(btype, bat); - if ((nr = get_atomtype_type(type, at)) != NOTSET) + if ((nr = at->atomTypeFromName(type)) != NOTSET) { auto message = gmx::formatString("Overriding atomtype %s", type); warning(wi, message); - if ((nr = set_atomtype(nr, at, symtab, atom, type, param, batype_nr, - atomnr)) == NOTSET) + if ((nr = at->setType(nr, symtab, *atom, type, param, batype_nr, + atomnr)) == NOTSET) { auto message = gmx::formatString("Replacing atomtype %s failed", type); warning_error_and_exit(wi, message, FARGS); } } - else if ((add_atomtype(at, symtab, atom, type, param, - batype_nr, atomnr)) == NOTSET) + else if ((at->addType(symtab, *atom, type, param, + batype_nr, atomnr)) == NOTSET) { auto message = gmx::formatString("Adding atomtype %s failed", type); warning_error_and_exit(wi, message, FARGS); @@ -715,10 +718,13 @@ static void push_bondtype(InteractionTypeParameters * bt, } } -void push_bt(Directive d, gmx::ArrayRef bt, int nral, - gpp_atomtype *at, - gpp_bond_atomtype *bat, char *line, - warninp *wi) +void push_bt(Directive d, + gmx::ArrayRef bt, + int nral, + PreprocessingAtomTypes *at, + gpp_bond_atomtype *bat, + char *line, + warninp *wi) { const char *formal[MAXATOMLIST+1] = { "%s", @@ -796,7 +802,7 @@ void push_bt(Directive d, gmx::ArrayRef bt, int nral, } for (i = 0; (i < nral); i++) { - if (at && ((p.a[i] = get_atomtype_type(alc[i], at)) == NOTSET)) + if ((at != nullptr) && ((p.a[i] = at->atomTypeFromName(alc[i])) == NOTSET)) { auto message = gmx::formatString("Unknown atomtype %s\n", alc[i]); warning_error_and_exit(wi, message, FARGS); @@ -984,7 +990,7 @@ void push_dihedraltype(Directive d, gmx::ArrayRef bt, } -void push_nbt(Directive d, t_nbparam **nbt, gpp_atomtype *atype, +void push_nbt(Directive d, t_nbparam **nbt, PreprocessingAtomTypes *atypes, char *pline, int nb_funct, warninp *wi) { @@ -1072,12 +1078,12 @@ void push_nbt(Directive d, t_nbparam **nbt, gpp_atomtype *atype, } /* Put the parameters in the matrix */ - if ((ai = get_atomtype_type (a0, atype)) == NOTSET) + if ((ai = atypes->atomTypeFromName(a0)) == NOTSET) { auto message = gmx::formatString("Atomtype %s not found", a0); warning_error_and_exit(wi, message, FARGS); } - if ((aj = get_atomtype_type (a1, atype)) == NOTSET) + if ((aj = atypes->atomTypeFromName(a1)) == NOTSET) { auto message = gmx::formatString("Atomtype %s not found", a1); warning_error_and_exit(wi, message, FARGS); @@ -1111,13 +1117,17 @@ void push_nbt(Directive d, t_nbparam **nbt, gpp_atomtype *atype, } void -push_cmaptype(Directive d, gmx::ArrayRef bt, int nral, gpp_atomtype *at, - gpp_bond_atomtype *bat, char *line, - warninp *wi) +push_cmaptype(Directive d, + gmx::ArrayRef bt, + int nral, + PreprocessingAtomTypes *at, + gpp_bond_atomtype *bat, + char *line, + warninp *wi) { const char *formal = "%s%s%s%s%s%s%s%s%n"; - int i, ft, ftype, nn, nrfp, nrfpA, nrfpB; + int ft, ftype, nn, nrfp, nrfpA, nrfpB; int start, nchar_consumed; int nxcmap, nycmap, ncmap, read_cmap, sl, nct; char s[20], alc[MAXATOMLIST+2][20]; @@ -1157,7 +1167,7 @@ push_cmaptype(Directive d, gmx::ArrayRef bt, int nral /* Read in CMAP parameters */ sl = 0; - for (i = 0; i < ncmap; i++) + for (int i = 0; i < ncmap; i++) { while (isspace(*(line+start+sl))) { @@ -1182,7 +1192,7 @@ push_cmaptype(Directive d, gmx::ArrayRef bt, int nral /* Check do that we got the number of parameters we expected */ if (read_cmap == nrfpA) { - for (i = 0; i < ncmap; i++) + for (int i = 0; i < ncmap; i++) { bt[F_CMAP].cmap.emplace_back(bt[F_CMAP].cmap[i]); } @@ -1211,7 +1221,7 @@ push_cmaptype(Directive d, gmx::ArrayRef bt, int nral bt[F_CMAP].cmapAngles++; /* Since we are incrementing here, we need to subtract later, see (*****) */ nct = (nral+1) * bt[F_CMAP].cmapAngles; - for (i = 0; (i < nral); i++) + for (int i = 0; (i < nral); i++) { if (at && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET)) { @@ -1239,7 +1249,7 @@ push_cmaptype(Directive d, gmx::ArrayRef bt, int nral } /* Is this correct?? */ - for (i = 0; i < MAXFORCEPARAM; i++) + for (int i = 0; i < MAXFORCEPARAM; i++) { p.c[i] = NOTSET; } @@ -1351,7 +1361,7 @@ static void push_cg(t_block *block, int *lastindex, int index, int a) } void push_atom(t_symtab *symtab, t_block *cgs, - t_atoms *at, gpp_atomtype *atype, char *line, int *lastcg, + t_atoms *at, PreprocessingAtomTypes *atypes, char *line, int *lastcg, warninp *wi) { int nr, ptype; @@ -1372,16 +1382,16 @@ void push_atom(t_symtab *symtab, t_block *cgs, return; } sscanf(id, "%d", &atomnr); - if ((type = get_atomtype_type(ctype, atype)) == NOTSET) + if ((type = atypes->atomTypeFromName(ctype)) == NOTSET) { auto message = gmx::formatString("Atomtype %s not found", ctype); warning_error_and_exit(wi, message, FARGS); } - ptype = get_atomtype_ptype(type, atype); + ptype = atypes->atomParticleTypeFromAtomType(type); /* Set default from type */ - q0 = get_atomtype_qA(type, atype); - m0 = get_atomtype_massA(type, atype); + q0 = atypes->atomChargeAFromAtomType(type); + m0 = atypes->atomMassAFromAtomType(type); typeB = type; qB = q0; mB = m0; @@ -1399,13 +1409,13 @@ void push_atom(t_symtab *symtab, t_block *cgs, m0 = mB = m; if (nscan > 2) { - if ((typeB = get_atomtype_type(ctypeB, atype)) == NOTSET) + if ((typeB = atypes->atomTypeFromName(ctypeB)) == NOTSET) { auto message = gmx::formatString("Atomtype %s not found", ctypeB); warning_error_and_exit(wi, message, FARGS); } - qB = get_atomtype_qA(typeB, atype); - mB = get_atomtype_massA(typeB, atype); + qB = atypes->atomChargeAFromAtomType(typeB); + mB = atypes->atomMassAFromAtomType(typeB); if (nscan > 3) { qB = qb; @@ -1424,7 +1434,7 @@ void push_atom(t_symtab *symtab, t_block *cgs, push_cg(cgs, lastcg, cgnumber, nr); - push_atom_now(symtab, at, atomnr, get_atomtype_atomnumber(type, atype), + push_atom_now(symtab, at, atomnr, atypes->atomNumberFromAtomType(type), type, ctype, ptype, resnumberic, resname, name, m0, q0, typeB, typeB == type ? ctype : ctypeB, mB, qB, wi); @@ -1558,20 +1568,20 @@ static bool default_nb_params(int ftype, gmx::ArrayRef bondtype, - t_atoms *at, gpp_atomtype *atype, + t_atoms *at, PreprocessingAtomTypes *atypes, t_param *p, bool bB, int *cmap_type, int *nparam_def, warninp *wi) { - int i, nparam_found; + int nparam_found; int ct; - bool bFound = FALSE; + bool bFound = false; nparam_found = 0; ct = 0; /* Match the current cmap angle against the list of cmap_types */ - for (i = 0; i < bondtype[F_CMAP].nct() && !bFound; i += 6) + for (int i = 0; i < bondtype[F_CMAP].nct() && !bFound; i += 6) { if (bB) { @@ -1580,14 +1590,14 @@ static bool default_cmap_params(gmx::ArrayRef bondtyp else { if ( - (get_atomtype_batype(at->atom[p->a[0]].type, atype) == bondtype[F_CMAP].cmapAtomTypes[i]) && - (get_atomtype_batype(at->atom[p->a[1]].type, atype) == bondtype[F_CMAP].cmapAtomTypes[i+1]) && - (get_atomtype_batype(at->atom[p->a[2]].type, atype) == bondtype[F_CMAP].cmapAtomTypes[i+2]) && - (get_atomtype_batype(at->atom[p->a[3]].type, atype) == bondtype[F_CMAP].cmapAtomTypes[i+3]) && - (get_atomtype_batype(at->atom[p->a[4]].type, atype) == bondtype[F_CMAP].cmapAtomTypes[i+4])) + (atypes->bondAtomTypeFromAtomType(at->atom[p->a[0]].type) == bondtype[F_CMAP].cmapAtomTypes[i]) && + (atypes->bondAtomTypeFromAtomType(at->atom[p->a[1]].type) == bondtype[F_CMAP].cmapAtomTypes[i+1]) && + (atypes->bondAtomTypeFromAtomType(at->atom[p->a[2]].type) == bondtype[F_CMAP].cmapAtomTypes[i+2]) && + (atypes->bondAtomTypeFromAtomType(at->atom[p->a[3]].type) == bondtype[F_CMAP].cmapAtomTypes[i+3]) && + (atypes->bondAtomTypeFromAtomType(at->atom[p->a[4]].type) == bondtype[F_CMAP].cmapAtomTypes[i+4])) { /* Found cmap torsion */ - bFound = TRUE; + bFound = true; ct = bondtype[F_CMAP].cmapAtomTypes[i+5]; nparam_found = 1; } @@ -1613,12 +1623,12 @@ static bool default_cmap_params(gmx::ArrayRef bondtyp */ static int natom_match(t_param *pi, int type_i, int type_j, int type_k, int type_l, - const gpp_atomtype* atype) + const PreprocessingAtomTypes* atypes) { - if ((pi->ai() == -1 || get_atomtype_batype(type_i, atype) == pi->ai()) && - (pi->aj() == -1 || get_atomtype_batype(type_j, atype) == pi->aj()) && - (pi->ak() == -1 || get_atomtype_batype(type_k, atype) == pi->ak()) && - (pi->al() == -1 || get_atomtype_batype(type_l, atype) == pi->al())) + if ((pi->ai() == -1 || atypes->bondAtomTypeFromAtomType(type_i) == pi->ai()) && + (pi->aj() == -1 || atypes->bondAtomTypeFromAtomType(type_j) == pi->aj()) && + (pi->ak() == -1 || atypes->bondAtomTypeFromAtomType(type_k) == pi->ak()) && + (pi->al() == -1 || atypes->bondAtomTypeFromAtomType(type_l) == pi->al())) { return (pi->ai() == -1 ? 0 : 1) + @@ -1633,7 +1643,7 @@ static int natom_match(t_param *pi, } static bool defaultInteractionTypeParameters(int ftype, gmx::ArrayRef bt, - t_atoms *at, gpp_atomtype *atype, + t_atoms *at, PreprocessingAtomTypes *atypes, t_param *p, bool bB, t_param **param_def, int *nparam_def) @@ -1672,11 +1682,11 @@ static bool defaultInteractionTypeParameters(int ftype, gmx::ArrayRefatom[p->ai()].typeB, at->atom[p->aj()].typeB, at->atom[p->ak()].typeB, at->atom[p->al()].typeB, atype); + nmatch = natom_match(pt, at->atom[p->ai()].typeB, at->atom[p->aj()].typeB, at->atom[p->ak()].typeB, at->atom[p->al()].typeB, atypes); } else { - nmatch = natom_match(pt, at->atom[p->ai()].type, at->atom[p->aj()].type, at->atom[p->ak()].type, at->atom[p->al()].type, atype); + nmatch = natom_match(pt, at->atom[p->ai()].type, at->atom[p->aj()].type, at->atom[p->ak()].type, at->atom[p->al()].type, atypes); } if (nmatch > nmatch_max) { @@ -1722,7 +1732,7 @@ static bool defaultInteractionTypeParameters(int ftype, gmx::ArrayRefatom[p->a[j]].typeB, atype) == pi->a[j])); j++) + (atypes->bondAtomTypeFromAtomType(at->atom[p->a[j]].typeB) == pi->a[j])); j++) { ; } @@ -1730,7 +1740,7 @@ static bool defaultInteractionTypeParameters(int ftype, gmx::ArrayRefatom[p->a[j]].type, atype) == pi->a[j])); j++) + (atypes->bondAtomTypeFromAtomType(at->atom[p->a[j]].type) == pi->a[j])); j++) { ; } @@ -1753,7 +1763,7 @@ static bool defaultInteractionTypeParameters(int ftype, gmx::ArrayRef bondtype, gmx::ArrayRef bond, - t_atoms *at, gpp_atomtype *atype, char *line, + t_atoms *at, PreprocessingAtomTypes *atypes, char *line, bool bBonded, bool bGenPairs, real fudgeQQ, bool bZero, bool *bWarn_copy_A_B, warninp *wi) @@ -1898,7 +1908,7 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, if (bBonded) { - bFoundA = defaultInteractionTypeParameters(ftype, bondtype, at, atype, ¶m, FALSE, ¶m_defA, &nparam_defA); + bFoundA = defaultInteractionTypeParameters(ftype, bondtype, at, atypes, ¶m, FALSE, ¶m_defA, &nparam_defA); if (bFoundA) { /* Copy the A-state and B-state default parameters. */ @@ -1908,7 +1918,7 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, param.c[j] = param_defA->c[j]; } } - bFoundB = defaultInteractionTypeParameters(ftype, bondtype, at, atype, ¶m, TRUE, ¶m_defB, &nparam_defB); + bFoundB = defaultInteractionTypeParameters(ftype, bondtype, at, atypes, ¶m, TRUE, ¶m_defB, &nparam_defB); if (bFoundB) { /* Copy only the B-state default parameters */ @@ -2162,7 +2172,7 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, void push_cmap(Directive d, gmx::ArrayRef bondtype, gmx::ArrayRef bond, - t_atoms *at, gpp_atomtype *atype, char *line, + t_atoms *at, PreprocessingAtomTypes *atypes, char *line, warninp *wi) { const char *aaformat[MAXATOMLIST+1] = @@ -2234,7 +2244,7 @@ void push_cmap(Directive d, gmx::ArrayRef bondtype, } /* Get the cmap type for this cmap angle */ - bFound = default_cmap_params(bondtype, at, atype, ¶m, FALSE, &cmap_type, &ncmap_params, wi); + bFound = default_cmap_params(bondtype, at, atypes, ¶m, FALSE, &cmap_type, &ncmap_params, wi); /* We want exactly one parameter (the cmap type in state A (currently no state B) back */ if (bFound && ncmap_params == 1) @@ -2476,7 +2486,7 @@ void push_excl(char *line, gmx::ArrayRef b2, warninp *wi) while (n == 1); } -int add_atomtype_decoupled(t_symtab *symtab, gpp_atomtype *at, +int add_atomtype_decoupled(t_symtab *symtab, PreprocessingAtomTypes *at, t_nbparam ***nbparam, t_nbparam ***pair) { t_atom atom; @@ -2495,7 +2505,7 @@ int add_atomtype_decoupled(t_symtab *symtab, gpp_atomtype *at, param.c[i] = 0.0; } - nr = add_atomtype(at, symtab, &atom, "decoupled", ¶m, -1, 0); + nr = at->addType(symtab, atom, "decoupled", ¶m, -1, 0); /* Add space in the non-bonded parameters matrix */ realloc_nb_params(at, nbparam, pair); diff --git a/src/gromacs/gmxpreprocess/toppush.h b/src/gromacs/gmxpreprocess/toppush.h index 48d706871a..33371cafc7 100644 --- a/src/gromacs/gmxpreprocess/toppush.h +++ b/src/gromacs/gmxpreprocess/toppush.h @@ -44,7 +44,7 @@ #include "gromacs/utility/real.h" enum class Directive : int; -struct gpp_atomtype; +class PreprocessingAtomTypes; struct gpp_bond_atomtype; struct t_atoms; struct t_block; @@ -61,41 +61,41 @@ struct ExclusionBlock; } // namespace gmx void generate_nbparams(int comb, int funct, InteractionTypeParameters *plist, - gpp_atomtype *atype, + PreprocessingAtomTypes *atype, warninp *wi); -void push_at (struct t_symtab *symtab, gpp_atomtype *at, +void push_at (struct t_symtab *symtab, PreprocessingAtomTypes *at, gpp_bond_atomtype *bat, char *line, int nb_funct, t_nbparam ***nbparam, t_nbparam ***pair, warninp *wi); void push_bt(Directive d, gmx::ArrayRef bt, int nral, - gpp_atomtype *at, gpp_bond_atomtype *bat, char *line, + PreprocessingAtomTypes *at, gpp_bond_atomtype *bat, char *line, warninp *wi); void push_dihedraltype(Directive d, gmx::ArrayRef bt, gpp_bond_atomtype *bat, char *line, warninp *wi); -void push_cmaptype(Directive d, gmx::ArrayRef bt, int nral, gpp_atomtype *at, +void push_cmaptype(Directive d, gmx::ArrayRef bt, int nral, PreprocessingAtomTypes *at, gpp_bond_atomtype *bat, char *line, warninp *wi); -void push_nbt(Directive d, t_nbparam **nbt, gpp_atomtype *atype, +void push_nbt(Directive d, t_nbparam **nbt, PreprocessingAtomTypes *atype, char *plines, int nb_funct, warninp *wi); -void push_atom(struct t_symtab *symtab, - t_block *cgs, - t_atoms *at, - gpp_atomtype *atype, - char *line, - int *lastcg, - warninp *wi); +void push_atom(struct t_symtab *symtab, + t_block *cgs, + t_atoms *at, + PreprocessingAtomTypes *atype, + char *line, + int *lastcg, + warninp *wi); void push_bond(Directive d, gmx::ArrayRef bondtype, gmx::ArrayRef bond, - t_atoms *at, gpp_atomtype *atype, char *line, + t_atoms *at, PreprocessingAtomTypes *atype, char *line, bool bBonded, bool bGenPairs, real fudgeQQ, bool bZero, bool *bWarn_copy_A_B, warninp *wi); @@ -103,7 +103,7 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, void push_cmap(Directive d, gmx::ArrayRef bondtype, gmx::ArrayRef bond, - t_atoms *at, gpp_atomtype *atype, char *line, + t_atoms *at, PreprocessingAtomTypes *atype, char *line, warninp *wi); void push_vsitesn(Directive d, gmx::ArrayRef bond, @@ -123,7 +123,7 @@ int copy_nbparams(t_nbparam **param, int ftype, InteractionTypeParameters *plist void free_nbparam(t_nbparam **param, int nr); -int add_atomtype_decoupled(struct t_symtab *symtab, gpp_atomtype *at, +int add_atomtype_decoupled(struct t_symtab *symtab, PreprocessingAtomTypes *at, t_nbparam ***nbparam, t_nbparam ***pair); /* Add an atom type with all parameters set to zero (no interactions). * Returns the atom type number. diff --git a/src/gromacs/gmxpreprocess/topshake.cpp b/src/gromacs/gmxpreprocess/topshake.cpp index 49dea59eee..ca0b56ff7b 100644 --- a/src/gromacs/gmxpreprocess/topshake.cpp +++ b/src/gromacs/gmxpreprocess/topshake.cpp @@ -83,7 +83,7 @@ static void copy_bond (InteractionTypeParameters *pr, int to, int from) static int count_hydrogens (char ***atomname, int nra, const int a[]) { - int i, nh; + int nh; if (!atomname) { @@ -92,7 +92,7 @@ static int count_hydrogens (char ***atomname, int nra, const int a[]) } nh = 0; - for (i = 0; (i < nra); i++) + for (int i = 0; (i < nra); i++) { if (toupper(**(atomname[a[i]])) == 'H') { @@ -105,13 +105,9 @@ static int count_hydrogens (char ***atomname, int nra, const int a[]) void make_shake(gmx::ArrayRef plist, t_atoms *atoms, int nshake) { - char ***info = atoms->atomname; - InteractionTypeParameters *pr; - InteractionTypeParameters *bonds; - t_param p, *bond, *ang; - real b_ij, b_jk; - int i, j, ftype, ftype_a; - bool bFound; + char ***info = atoms->atomname; + real b_ij, b_jk; + int i, j; if (nshake != eshNONE) { @@ -138,27 +134,25 @@ void make_shake(gmx::ArrayRef plist, t_atoms *atoms, /* Add all the angles with hydrogens to the shake list * and remove them from the bond list */ - for (ftype = 0; (ftype < F_NRE); ftype++) + for (int ftype = 0; (ftype < F_NRE); ftype++) { if (interaction_function[ftype].flags & IF_BTYPE) { - bonds = &(plist[ftype]); + InteractionTypeParameters *bonds = &(plist[ftype]); - for (ftype_a = 0; (bonds->nr > 0 && ftype_a < F_NRE); ftype_a++) + for (int ftype_a = 0; (bonds->nr > 0 && ftype_a < F_NRE); ftype_a++) { if (interaction_function[ftype_a].flags & IF_ATYPE) { - pr = &(plist[ftype_a]); + InteractionTypeParameters *pr = &(plist[ftype_a]); - for (i = 0; (i < pr->nr); ) + for (int i = 0; (i < pr->nr); ) { - int numhydrogens; - - ang = &(pr->param[i]); + t_param *ang = &(pr->param[i]); #ifdef DEBUG printf("Angle: %d-%d-%d\n", ang->ai(), ang->aj(), ang->ak()); #endif - numhydrogens = count_hydrogens(info, 3, ang->a); + int numhydrogens = count_hydrogens(info, 3, ang->a); if ((nshake == eshALLANGLES) || (numhydrogens > 1) || (numhydrogens == 1 && toupper(**(info[ang->a[1]])) == 'O')) @@ -167,15 +161,16 @@ void make_shake(gmx::ArrayRef plist, t_atoms *atoms, * are constrained. * append this angle to the shake list */ + t_param p; p.ai() = ang->ai(); p.aj() = ang->ak(); /* Calculate length of constraint */ - bFound = FALSE; + bool bFound = false; b_ij = b_jk = 0.0; - for (j = 0; !bFound && (j < bonds->nr); j++) + for (int j = 0; !bFound && (j < bonds->nr); j++) { - bond = &(bonds->param[j]); + t_param *bond = &(bonds->param[j]); if (((bond->ai() == ang->ai()) && (bond->aj() == ang->aj())) || ((bond->ai() == ang->aj()) && @@ -222,11 +217,11 @@ void make_shake(gmx::ArrayRef plist, t_atoms *atoms, /* Add all the bonds with hydrogens to the shake list * and remove them from the bond list */ - for (ftype = 0; (ftype < F_NRE); ftype++) + for (int ftype = 0; (ftype < F_NRE); ftype++) { if (interaction_function[ftype].flags & IF_BTYPE) { - pr = &(plist[ftype]); + InteractionTypeParameters *pr = &(plist[ftype]); j = 0; for (i = 0; i < pr->nr; i++) { @@ -234,6 +229,7 @@ void make_shake(gmx::ArrayRef plist, t_atoms *atoms, (count_hydrogens (info, 2, pr->param[i].a) > 0) ) { /* append this bond to the shake list */ + t_param p; p.ai() = pr->param[i].ai(); p.aj() = pr->param[i].aj(); p.c0() = pr->param[i].c0(); diff --git a/src/gromacs/gmxpreprocess/toputil.cpp b/src/gromacs/gmxpreprocess/toputil.cpp index ecde0d9fa8..5fc3e04b0f 100644 --- a/src/gromacs/gmxpreprocess/toputil.cpp +++ b/src/gromacs/gmxpreprocess/toputil.cpp @@ -131,7 +131,7 @@ void done_plist(gmx::ArrayRef plist) } } -void cp_param(t_param *dest, t_param *src) +void cp_param(t_param *dest, const t_param *src) { int j; @@ -148,17 +148,15 @@ void cp_param(t_param *dest, t_param *src) void add_param_to_list(InteractionTypeParameters *list, t_param *b) { - int j; - /* allocate one position extra */ pr_alloc (1, list); /* fill the arrays */ - for (j = 0; (j < MAXFORCEPARAM); j++) + for (int j = 0; (j < MAXFORCEPARAM); j++) { list->param[list->nr].c[j] = b->c[j]; } - for (j = 0; (j < MAXATOMLIST); j++) + for (int j = 0; (j < MAXATOMLIST); j++) { list->param[list->nr].a[j] = b->a[j]; } @@ -169,7 +167,7 @@ void add_param_to_list(InteractionTypeParameters *list, t_param *b) /* PRINTING STRUCTURES */ -static void print_bt(FILE *out, Directive d, gpp_atomtype *at, +static void print_bt(FILE *out, Directive d, PreprocessingAtomTypes *at, int ftype, int fsubtype, gmx::ArrayRef plist, bool bFullDih) { @@ -177,17 +175,17 @@ static void print_bt(FILE *out, Directive d, gpp_atomtype *at, * all four atoms to determine the type. */ const int dihp[2][2] = { { 1, 2 }, { 0, 3 } }; - int i, j, f, nral, nrfp; - bool bDih = FALSE, bSwapParity; + int nral, nrfp; + bool bDih = false, bSwapParity; const InteractionTypeParameters *bt = &(plist[ftype]); - if (!bt->nr) + if (bt->nr == 0) { return; } - f = 0; + int f = 0; switch (ftype) { case F_G96ANGLES: @@ -266,42 +264,42 @@ static void print_bt(FILE *out, Directive d, gpp_atomtype *at, if (!bDih) { fprintf (out, "%3s %4s", "ai", "aj"); - for (j = 2; (j < nral); j++) + for (int j = 2; (j < nral); j++) { fprintf (out, " %3c%c", 'a', 'i'+j); } } else { - for (j = 0; (j < 2); j++) + for (int j = 0; (j < 2); j++) { fprintf (out, "%3c%c", 'a', 'i'+dihp[f][j]); } } fprintf (out, " funct"); - for (j = 0; (j < nrfp); j++) + for (int j = 0; (j < nrfp); j++) { fprintf (out, " %12c%1d", 'c', j); } fprintf (out, "\n"); /* print bondtypes */ - for (i = 0; (i < bt->nr); i++) + for (int i = 0; (i < bt->nr); i++) { bSwapParity = (bt->param[i].c0() == NOTSET) && (bt->param[i].c1() == -1); if (!bDih) { - for (j = 0; (j < nral); j++) + for (int j = 0; (j < nral); j++) { - fprintf (out, "%5s ", get_atomtype_name(bt->param[i].a[j], at)); + fprintf (out, "%5s ", at->atomNameFromAtomType(bt->param[i].a[j])); } } else { - for (j = 0; (j < 2); j++) + for (int j = 0; (j < 2); j++) { - fprintf (out, "%5s ", get_atomtype_name(bt->param[i].a[dihp[f][j]], at)); + fprintf (out, "%5s ", at->atomNameFromAtomType(bt->param[i].a[dihp[f][j]])); } } fprintf (out, "%5d ", bSwapParity ? -f-1 : f+1); @@ -312,7 +310,7 @@ static void print_bt(FILE *out, Directive d, gpp_atomtype *at, } else { - for (j = 0; (j < nrfp && (bt->param[i].c[j] != NOTSET)); j++) + for (int j = 0; (j < nrfp && (bt->param[i].c[j] != NOTSET)); j++) { fprintf (out, "%13.6e ", bt->param[i].c[j]); } @@ -396,14 +394,14 @@ static double get_residue_charge(const t_atoms *atoms, int at) return q; } -void print_atoms(FILE *out, gpp_atomtype *atype, t_atoms *at, int *cgnr, +void print_atoms(FILE *out, PreprocessingAtomTypes *atype, t_atoms *at, int *cgnr, bool bRTPresname) { - int i, ri; - int tpA, tpB; - const char *as; - char *tpnmA, *tpnmB; - double qres, qtot; + int i, ri; + int tpA, tpB; + const char *as; + const char *tpnmA, *tpnmB; + double qres, qtot; as = dir2str(Directive::d_atoms); fprintf(out, "[ %s ]\n", as); @@ -437,7 +435,7 @@ void print_atoms(FILE *out, gpp_atomtype *atype, t_atoms *at, int *cgnr, fprintf(out, "\n"); } tpA = at->atom[i].type; - if ((tpnmA = get_atomtype_name(tpA, atype)) == nullptr) + if ((tpnmA = atype->atomNameFromAtomType(tpA)) == nullptr) { gmx_fatal(FARGS, "tpA = %d, i= %d in print_atoms", tpA, i); } @@ -456,7 +454,7 @@ void print_atoms(FILE *out, gpp_atomtype *atype, t_atoms *at, int *cgnr, if (PERTURBED(at->atom[i])) { tpB = at->atom[i].typeB; - if ((tpnmB = get_atomtype_name(tpB, atype)) == nullptr) + if ((tpnmB = atype->atomNameFromAtomType(tpB)) == nullptr) { gmx_fatal(FARGS, "tpB = %d, i= %d in print_atoms", tpB, i); } @@ -491,26 +489,23 @@ void print_atoms(FILE *out, gpp_atomtype *atype, t_atoms *at, int *cgnr, void print_bondeds(FILE *out, int natoms, Directive d, int ftype, int fsubtype, gmx::ArrayRef plist) { - t_symtab stab; - gpp_atomtype *atype; - t_param *param; - t_atom *a; - int i; + t_symtab stab; + t_param *param; + t_atom *a; - atype = init_atomtype(); + PreprocessingAtomTypes atype; snew(a, 1); snew(param, 1); open_symtab(&stab); - for (i = 0; (i < natoms); i++) + for (int i = 0; (i < natoms); i++) { char buf[12]; sprintf(buf, "%4d", (i+1)); - add_atomtype(atype, &stab, a, buf, param, 0, 0); + atype.addType(&stab, *a, buf, param, 0, 0); } - print_bt(out, d, atype, ftype, fsubtype, plist, TRUE); + print_bt(out, d, &atype, ftype, fsubtype, plist, TRUE); done_symtab(&stab); sfree(a); sfree(param); - done_atomtype(atype); } diff --git a/src/gromacs/gmxpreprocess/toputil.h b/src/gromacs/gmxpreprocess/toputil.h index 068225238e..2891177b17 100644 --- a/src/gromacs/gmxpreprocess/toputil.h +++ b/src/gromacs/gmxpreprocess/toputil.h @@ -42,8 +42,10 @@ #include "gromacs/utility/arrayref.h" +#include "gromacs/utility/arrayref.h" + enum class Directive : int; -struct gpp_atomtype; +class PreprocessingAtomTypes; struct t_atoms; struct t_blocka; struct t_excls; @@ -59,7 +61,7 @@ void pr_alloc (int extra, InteractionTypeParameters *pr); void set_p_string(t_param *p, const char *s); -void cp_param(t_param *dest, t_param *src); +void cp_param(t_param *dest, const t_param *src); void add_param_to_list(InteractionTypeParameters *list, t_param *b); @@ -74,7 +76,7 @@ void done_plist(gmx::ArrayRef plist); void print_blocka(FILE *out, const char *szName, const char *szIndex, const char *szA, t_blocka *block); -void print_atoms(FILE *out, gpp_atomtype *atype, t_atoms *at, int *cgnr, +void print_atoms(FILE *out, PreprocessingAtomTypes *atype, t_atoms *at, int *cgnr, bool bRTPresname); void print_bondeds(FILE *out, int natoms, Directive d, diff --git a/src/gromacs/gmxpreprocess/vsite_parm.cpp b/src/gromacs/gmxpreprocess/vsite_parm.cpp index a20d19353d..c0c3b36dbd 100644 --- a/src/gromacs/gmxpreprocess/vsite_parm.cpp +++ b/src/gromacs/gmxpreprocess/vsite_parm.cpp @@ -372,11 +372,9 @@ static real get_angle(int nrang, t_mybonded angles[], return angle; } -static char *get_atomtype_name_AB(t_atom *atom, gpp_atomtype *atype) +static const char *get_atomtype_name_AB(t_atom *atom, PreprocessingAtomTypes *atypes) { - char *name; - - name = get_atomtype_name(atom->type, atype); + const char* name = atypes->atomNameFromAtomType(atom->type); /* When using the decoupling option, atom types are changed * to decoupled for the non-bonded interactions, but the virtual @@ -390,13 +388,13 @@ static char *get_atomtype_name_AB(t_atom *atom, gpp_atomtype *atype) */ if (strcmp(name, "decoupled") == 0) { - name = get_atomtype_name(atom->typeB, atype); + name = atypes->atomNameFromAtomType(atom->typeB); } return name; } -static bool calc_vsite3_param(gpp_atomtype *atype, +static bool calc_vsite3_param(PreprocessingAtomTypes *atypes, t_param *param, t_atoms *at, int nrbond, t_mybonded *bonds, int nrang, t_mybonded *angles ) @@ -411,10 +409,10 @@ static bool calc_vsite3_param(gpp_atomtype *atype, /* check if this is part of a NH3 , NH2-umbrella or CH3 group, * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */ bXH3 = - ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atype), "MNH", 3) == 0) && - (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atype), "MNH", 3) == 0) ) || - ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atype), "MCH3", 4) == 0) && - (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atype), "MCH3", 4) == 0) ); + ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atypes), "MNH", 3) == 0) && + (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atypes), "MNH", 3) == 0) ) || + ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atypes), "MCH3", 4) == 0) && + (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atypes), "MCH3", 4) == 0) ); bjk = get_bond_length(nrbond, bonds, param->aj(), param->ak()); bjl = get_bond_length(nrbond, bonds, param->aj(), param->al()); @@ -534,7 +532,7 @@ static bool calc_vsite3fad_param(t_param *param, return bError; } -static bool calc_vsite3out_param(gpp_atomtype *atype, +static bool calc_vsite3out_param(PreprocessingAtomTypes *atypes, t_param *param, t_atoms *at, int nrbond, t_mybonded *bonds, int nrang, t_mybonded *angles) @@ -551,10 +549,10 @@ static bool calc_vsite3out_param(gpp_atomtype *atype, /* check if this is part of a NH2-umbrella, NH3 or CH3 group, * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */ bXH3 = - ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atype), "MNH", 3) == 0) && - (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atype), "MNH", 3) == 0) ) || - ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atype), "MCH3", 4) == 0) && - (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atype), "MCH3", 4) == 0) ); + ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atypes), "MNH", 3) == 0) && + (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atypes), "MNH", 3) == 0) ) || + ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->ak()], atypes), "MCH3", 4) == 0) && + (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->al()], atypes), "MCH3", 4) == 0) ); /* check if construction parity must be swapped */ bSwapParity = ( param->c1() == -1 ); @@ -746,7 +744,7 @@ calc_vsite4fdn_param(t_param *param, -int set_vsites(bool bVerbose, t_atoms *atoms, gpp_atomtype *atype, +int set_vsites(bool bVerbose, t_atoms *atoms, PreprocessingAtomTypes *atypes, gmx::ArrayRef plist) { int i, j, ftype; @@ -817,7 +815,7 @@ int set_vsites(bool bVerbose, t_atoms *atoms, gpp_atomtype *atype, { case F_VSITE3: bERROR = - calc_vsite3_param(atype, &(plist[ftype].param[i]), atoms, + calc_vsite3_param(atypes, &(plist[ftype].param[i]), atoms, nrbond, bonds, nrang, angles); break; case F_VSITE3FD: @@ -832,7 +830,7 @@ int set_vsites(bool bVerbose, t_atoms *atoms, gpp_atomtype *atype, break; case F_VSITE3OUT: bERROR = - calc_vsite3out_param(atype, &(plist[ftype].param[i]), atoms, + calc_vsite3out_param(atypes, &(plist[ftype].param[i]), atoms, nrbond, bonds, nrang, angles); break; case F_VSITE4FD: diff --git a/src/gromacs/gmxpreprocess/vsite_parm.h b/src/gromacs/gmxpreprocess/vsite_parm.h index 81005a0fe4..ec383ac64c 100644 --- a/src/gromacs/gmxpreprocess/vsite_parm.h +++ b/src/gromacs/gmxpreprocess/vsite_parm.h @@ -40,12 +40,12 @@ #include "gromacs/utility/arrayref.h" -struct gpp_atomtype; +class PreprocessingAtomTypes; struct gmx_moltype_t; struct t_atoms; struct InteractionTypeParameters; -int set_vsites(bool bVerbose, t_atoms *atoms, gpp_atomtype *atype, +int set_vsites(bool bVerbose, t_atoms *atoms, PreprocessingAtomTypes *atype, gmx::ArrayRef plist); /* set parameters for virtual sites, return number of virtual sites */ diff --git a/src/gromacs/gmxpreprocess/x2top.cpp b/src/gromacs/gmxpreprocess/x2top.cpp index a120b42925..3a587ab5c3 100644 --- a/src/gromacs/gmxpreprocess/x2top.cpp +++ b/src/gromacs/gmxpreprocess/x2top.cpp @@ -178,25 +178,26 @@ static int *set_cgnr(t_atoms *atoms, bool bUsePDBcharge, real *qtot, real *mtot) return cgnr; } -static gpp_atomtype *set_atom_type(t_symtab *tab, t_atoms *atoms, InteractionTypeParameters *bonds, - int *nbonds, int nnm, t_nm2type nm2t[]) +static void set_atom_type(PreprocessingAtomTypes *atypes, + t_symtab *tab, + t_atoms *atoms, + InteractionTypeParameters *bonds, + int *nbonds, + int nnm, + t_nm2type nm2t[]) { - gpp_atomtype *atype; int nresolved; - atype = init_atomtype(); snew(atoms->atomtype, atoms->nr); - nresolved = nm2type(nnm, nm2t, tab, atoms, atype, nbonds, bonds); + nresolved = nm2type(nnm, nm2t, tab, atoms, atypes, nbonds, bonds); if (nresolved != atoms->nr) { gmx_fatal(FARGS, "Could only find a forcefield type for %d out of %d atoms", nresolved, atoms->nr); } - fprintf(stderr, "There are %d different atom types in your sample\n", - get_atomtype_ntypes(atype)); - - return atype; + fprintf(stderr, "There are %zu different atom types in your sample\n", + atypes->size()); } static void lo_set_force_const(InteractionTypeParameters *plist, real c[], int nrfp, bool bRound, @@ -298,9 +299,7 @@ static void calc_angles_dihs(InteractionTypeParameters *ang, InteractionTypePara static void dump_hybridization(FILE *fp, t_atoms *atoms, int nbonds[]) { - int i; - - for (i = 0; (i < atoms->nr); i++) + for (int i = 0; (i < atoms->nr); i++) { fprintf(fp, "Atom %5s has %1d bonds\n", *atoms->atomname[i], nbonds[i]); } @@ -333,12 +332,16 @@ static void print_pl(FILE *fp, gmx::ArrayRef pl } } -static void print_rtp(const char *filenm, const char *title, t_atoms *atoms, - gmx::ArrayRef plist, gpp_atomtype *atype, int cgnr[]) +static void print_rtp(const char *filenm, + const char *title, + t_atoms *atoms, + gmx::ArrayRef plist, + PreprocessingAtomTypes *atypes, + int cgnr[]) { - FILE *fp; - int i, tp; - char *tpnm; + FILE *fp; + int i, tp; + const char *tpnm; fp = gmx_fio_fopen(filenm, "w"); fprintf(fp, "; %s\n", title); @@ -349,7 +352,7 @@ static void print_rtp(const char *filenm, const char *title, t_atoms *atoms, for (i = 0; (i < atoms->nr); i++) { tp = atoms->atom[i].type; - if ((tpnm = get_atomtype_name(tp, atype)) == nullptr) + if ((tpnm = atypes->atomNameFromAtomType(tp)) == nullptr) { gmx_fatal(FARGS, "tp = %d, i = %d in print_rtp", tp, i); } @@ -397,7 +400,6 @@ int gmx_x2top(int argc, char *argv[]) FILE *fp; std::array plist; t_excls *excls; - gpp_atomtype *atype; t_nextnb nnb; t_nm2type *nm2t; t_mols mymol; @@ -526,7 +528,8 @@ int gmx_x2top(int argc, char *argv[]) mk_bonds(nnm, nm2t, atoms, x, &(plist[F_BONDS]), nbonds, bPBC, box); open_symtab(&symtab); - atype = set_atom_type(&symtab, atoms, &(plist[F_BONDS]), nbonds, nnm, nm2t); + PreprocessingAtomTypes atypes; + set_atom_type(&atypes, &symtab, atoms, &(plist[F_BONDS]), nbonds, nnm, nm2t); /* Make Angles and Dihedrals */ snew(excls, atoms->nr); @@ -566,7 +569,7 @@ int gmx_x2top(int argc, char *argv[]) fp = ftp2FILE(efTOP, NFILE, fnm, "w"); print_top_header(fp, ftp2fn(efTOP, NFILE, fnm), TRUE, ffdir, 1.0); - write_top(fp, nullptr, mymol.name.c_str(), atoms, FALSE, bts, plist, excls, atype, + write_top(fp, nullptr, mymol.name.c_str(), atoms, FALSE, bts, plist, excls, &atypes, cgnr, rtp_header_settings.nrexcl); print_top_mols(fp, mymol.name.c_str(), ffdir, nullptr, {}, gmx::arrayRefFromArray(&mymol, 1)); @@ -575,7 +578,7 @@ int gmx_x2top(int argc, char *argv[]) if (bRTP) { print_rtp(ftp2fn(efRTP, NFILE, fnm), "Generated by x2top", - atoms, plist, atype, cgnr); + atoms, plist, &atypes, cgnr); } if (debug) -- 2.11.4.GIT