From 0f8d4458804477355483e25678697cc8fc05ddfb Mon Sep 17 00:00:00 2001 From: Paul Bauer Date: Fri, 22 Feb 2019 18:07:39 +0100 Subject: [PATCH] Refactor t_param to InteractionType Refs #2833 Change-Id: Iab043e96399c2a9615e66c5889010331c95782df --- src/gromacs/gmxpreprocess/add_par.cpp | 152 +--- src/gromacs/gmxpreprocess/add_par.h | 7 +- src/gromacs/gmxpreprocess/convparm.cpp | 27 +- src/gromacs/gmxpreprocess/gen_ad.cpp | 597 +++++--------- src/gromacs/gmxpreprocess/gen_vsite.cpp | 83 +- src/gromacs/gmxpreprocess/gpp_atomtype.cpp | 124 ++- src/gromacs/gmxpreprocess/gpp_atomtype.h | 28 +- src/gromacs/gmxpreprocess/gpp_nextnb.cpp | 10 +- src/gromacs/gmxpreprocess/grompp.cpp | 174 +++- src/gromacs/gmxpreprocess/grompp_impl.h | 126 ++- src/gromacs/gmxpreprocess/nm2type.cpp | 19 +- src/gromacs/gmxpreprocess/pdb2top.cpp | 107 +-- src/gromacs/gmxpreprocess/resall.cpp | 6 +- src/gromacs/gmxpreprocess/tomorse.cpp | 61 +- src/gromacs/gmxpreprocess/topio.cpp | 32 +- src/gromacs/gmxpreprocess/toppush.cpp | 889 +++++++++++---------- src/gromacs/gmxpreprocess/toppush.h | 2 +- src/gromacs/gmxpreprocess/topshake.cpp | 105 +-- src/gromacs/gmxpreprocess/toputil.cpp | 132 +-- src/gromacs/gmxpreprocess/toputil.h | 14 +- src/gromacs/gmxpreprocess/vsite_parm.cpp | 488 ++++++----- src/gromacs/gmxpreprocess/x2top.cpp | 88 +- 22 files changed, 1468 insertions(+), 1803 deletions(-) diff --git a/src/gromacs/gmxpreprocess/add_par.cpp b/src/gromacs/gmxpreprocess/add_par.cpp index a3b8283c65..046eecceb5 100644 --- a/src/gromacs/gmxpreprocess/add_par.cpp +++ b/src/gromacs/gmxpreprocess/add_par.cpp @@ -52,170 +52,80 @@ #include "hackblock.h" -static void clear_atom_list(int i0, int a[]) +void add_param(InteractionTypeParameters *ps, + int ai, + int aj, + gmx::ArrayRef c, + const char *s) { - int i; - - for (i = i0; i < MAXATOMLIST; i++) - { - a[i] = -1; - } -} - -static void clear_force_param(int i0, real c[]) -{ - int i; - - for (i = i0; i < MAXFORCEPARAM; i++) - { - c[i] = NOTSET; - } -} - -void add_param(InteractionTypeParameters *ps, int ai, int aj, const real *c, const char *s) -{ - int i; - if ((ai < 0) || (aj < 0)) { gmx_fatal(FARGS, "Trying to add impossible atoms: ai=%d, aj=%d", ai, aj); } - pr_alloc(1, ps); - ps->param[ps->nr].ai() = ai; - ps->param[ps->nr].aj() = aj; - clear_atom_list(2, ps->param[ps->nr].a); - if (c == nullptr) - { - clear_force_param(0, ps->param[ps->nr].c); - } - else - { - for (i = 0; (i < MAXFORCEPARAM); i++) - { - ps->param[ps->nr].c[i] = c[i]; - } - } - set_p_string(&(ps->param[ps->nr]), s); - ps->nr++; + std::vector atoms = {ai, aj}; + std::vector forceParm(c.begin(), c.end()); + + ps->interactionTypes.emplace_back(InteractionType(atoms, forceParm, s ? s : "")); } void add_imp_param(InteractionTypeParameters *ps, int ai, int aj, int ak, int al, real c0, real c1, const char *s) { - pr_alloc(1, ps); - ps->param[ps->nr].ai() = ai; - ps->param[ps->nr].aj() = aj; - ps->param[ps->nr].ak() = ak; - ps->param[ps->nr].al() = al; - clear_atom_list (4, ps->param[ps->nr].a); - ps->param[ps->nr].c0() = c0; - ps->param[ps->nr].c1() = c1; - clear_force_param(2, ps->param[ps->nr].c); - set_p_string(&(ps->param[ps->nr]), s); - ps->nr++; + std::vector atoms = {ai, aj, ak, al}; + std::vector forceParm = {c0, c1}; + ps->interactionTypes.emplace_back(InteractionType(atoms, forceParm, s ? s : "")); } void add_dih_param(InteractionTypeParameters *ps, int ai, int aj, int ak, int al, real c0, real c1, real c2, const char *s) { - pr_alloc(1, ps); - ps->param[ps->nr].ai() = ai; - ps->param[ps->nr].aj() = aj; - ps->param[ps->nr].ak() = ak; - ps->param[ps->nr].al() = al; - clear_atom_list (4, ps->param[ps->nr].a); - ps->param[ps->nr].c0() = c0; - ps->param[ps->nr].c1() = c1; - ps->param[ps->nr].c2() = c2; - clear_force_param(3, ps->param[ps->nr].c); - set_p_string(&(ps->param[ps->nr]), s); - ps->nr++; + std::vector atoms = {ai, aj, ak, al}; + std::vector forceParm = {c0, c1, c2}; + ps->interactionTypes.emplace_back(InteractionType(atoms, forceParm, s ? s : "")); } void add_cmap_param(InteractionTypeParameters *ps, int ai, int aj, int ak, int al, int am, const char *s) { - pr_alloc(1, ps); - ps->param[ps->nr].ai() = ai; - ps->param[ps->nr].aj() = aj; - ps->param[ps->nr].ak() = ak; - ps->param[ps->nr].al() = al; - ps->param[ps->nr].am() = am; - clear_atom_list(5, ps->param[ps->nr].a); - clear_force_param(0, ps->param[ps->nr].c); - set_p_string(&(ps->param[ps->nr]), s); - ps->nr++; + std::vector atoms = {ai, aj, ak, al, am}; + ps->interactionTypes.emplace_back(InteractionType(atoms, {}, s ? s : "")); } void add_vsite2_atoms(InteractionTypeParameters *ps, int ai, int aj, int ak) { - pr_alloc(1, ps); - ps->param[ps->nr].ai() = ai; - ps->param[ps->nr].aj() = aj; - ps->param[ps->nr].ak() = ak; - clear_atom_list (3, ps->param[ps->nr].a); - clear_force_param(0, ps->param[ps->nr].c); - set_p_string(&(ps->param[ps->nr]), ""); - ps->nr++; + std::vector atoms = {ai, aj, ak}; + ps->interactionTypes.emplace_back(InteractionType(atoms, {})); } void add_vsite2_param(InteractionTypeParameters *ps, int ai, int aj, int ak, real c0) { - pr_alloc(1, ps); - ps->param[ps->nr].ai() = ai; - ps->param[ps->nr].aj() = aj; - ps->param[ps->nr].ak() = ak; - clear_atom_list (3, ps->param[ps->nr].a); - ps->param[ps->nr].c0() = c0; - clear_force_param(1, ps->param[ps->nr].c); - set_p_string(&(ps->param[ps->nr]), ""); - ps->nr++; + std::vector atoms = {ai, aj, ak}; + std::vector forceParm = {c0}; + ps->interactionTypes.emplace_back(InteractionType(atoms, forceParm)); } void add_vsite3_param(InteractionTypeParameters *ps, int ai, int aj, int ak, int al, real c0, real c1) { - pr_alloc(1, ps); - ps->param[ps->nr].ai() = ai; - ps->param[ps->nr].aj() = aj; - ps->param[ps->nr].ak() = ak; - ps->param[ps->nr].al() = al; - clear_atom_list (4, ps->param[ps->nr].a); - ps->param[ps->nr].c0() = c0; - ps->param[ps->nr].c1() = c1; - clear_force_param(2, ps->param[ps->nr].c); - set_p_string(&(ps->param[ps->nr]), ""); - ps->nr++; + std::vector atoms = {ai, aj, ak, al}; + std::vector forceParm = {c0, c1}; + ps->interactionTypes.emplace_back(InteractionType(atoms, forceParm)); } void add_vsite3_atoms(InteractionTypeParameters *ps, int ai, int aj, int ak, int al, bool bSwapParity) { - pr_alloc(1, ps); - ps->param[ps->nr].ai() = ai; - ps->param[ps->nr].aj() = aj; - ps->param[ps->nr].ak() = ak; - ps->param[ps->nr].al() = al; - clear_atom_list (4, ps->param[ps->nr].a); - clear_force_param(0, ps->param[ps->nr].c); + std::vector atoms = {ai, aj, ak, al}; + ps->interactionTypes.emplace_back(InteractionType(atoms, {})); + if (bSwapParity) { - ps->param[ps->nr].c1() = -1; + ps->interactionTypes.back().setForceParameter(1, -1); } - set_p_string(&(ps->param[ps->nr]), ""); - ps->nr++; } void add_vsite4_atoms(InteractionTypeParameters *ps, int ai, int aj, int ak, int al, int am) { - pr_alloc(1, ps); - ps->param[ps->nr].ai() = ai; - ps->param[ps->nr].aj() = aj; - ps->param[ps->nr].ak() = ak; - ps->param[ps->nr].al() = al; - ps->param[ps->nr].am() = am; - clear_atom_list (5, ps->param[ps->nr].a); - clear_force_param(0, ps->param[ps->nr].c); - set_p_string(&(ps->param[ps->nr]), ""); - ps->nr++; + std::vector atoms = {ai, aj, ak, al, am}; + ps->interactionTypes.emplace_back(InteractionType(atoms, {})); } int search_jtype(const PreprocessResidue &localPpResidue, const char *name, bool bNterm) diff --git a/src/gromacs/gmxpreprocess/add_par.h b/src/gromacs/gmxpreprocess/add_par.h index de35358e4b..afc517c1be 100644 --- a/src/gromacs/gmxpreprocess/add_par.h +++ b/src/gromacs/gmxpreprocess/add_par.h @@ -38,12 +38,17 @@ #ifndef GMX_GMXPREPROCESS_ADD_PAR_H #define GMX_GMXPREPROCESS_ADD_PAR_H +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/real.h" struct InteractionTypeParameters; struct PreprocessResidue; -void add_param(InteractionTypeParameters *ps, int ai, int aj, const real *c, const char *s); +void add_param(InteractionTypeParameters *ps, + int ai, + int aj, + gmx::ArrayRef c, + const char *s); void add_imp_param(InteractionTypeParameters *ps, int ai, int aj, int ak, int al, real c0, real c1, const char *s); diff --git a/src/gromacs/gmxpreprocess/convparm.cpp b/src/gromacs/gmxpreprocess/convparm.cpp index 546cd97ae0..f3e6fb9c56 100644 --- a/src/gromacs/gmxpreprocess/convparm.cpp +++ b/src/gromacs/gmxpreprocess/convparm.cpp @@ -107,7 +107,7 @@ static void set_ljparams(int comb, double reppow, double v, double w, */ static int assign_param(t_functype ftype, t_iparams *newparam, - real old[MAXFORCEPARAM], int comb, double reppow) + gmx::ArrayRef old, int comb, double reppow) { bool all_param_zero = true; @@ -446,7 +446,7 @@ assign_param(t_functype ftype, t_iparams *newparam, } static int enter_params(gmx_ffparams_t *ffparams, t_functype ftype, - real forceparams[MAXFORCEPARAM], int comb, real reppow, + gmx::ArrayRef forceparams, int comb, real reppow, int start, bool bAppend) { t_iparams newparam; @@ -487,12 +487,12 @@ static int enter_params(gmx_ffparams_t *ffparams, t_functype ftype, } static void append_interaction(InteractionList *ilist, - int type, int nral, const int a[MAXATOMLIST]) + int type, gmx::ArrayRef a) { ilist->iatoms.push_back(type); - for (int i = 0; (i < nral); i++) + for (const auto &atom : a) { - ilist->iatoms.push_back(a[i]); + ilist->iatoms.push_back(atom); } } @@ -500,20 +500,17 @@ static void enter_function(const InteractionTypeParameters *p, t_functype ftype, gmx_ffparams_t *ffparams, InteractionList *il, bool bNB, bool bAppend) { - int k, type, nr, nral, start; + int start = ffparams->numTypes(); - start = ffparams->numTypes(); - nr = p->nr; - - for (k = 0; k < nr; k++) + for (auto &parm : p->interactionTypes) { - type = enter_params(ffparams, ftype, p->param[k].c, comb, reppow, start, bAppend); + int type = enter_params(ffparams, ftype, parm.forceParam(), comb, reppow, start, bAppend); /* Type==-1 is used as a signal that this interaction is all-zero and should not be added. */ if (!bNB && type >= 0) { - assert(il); - nral = NRAL(ftype); - append_interaction(il, type, nral, p->param[k].a); + GMX_RELEASE_ASSERT(il, "Need valid interaction list"); + GMX_RELEASE_ASSERT(parm.atoms().ssize() == NRAL(ftype), "Need to have correct number of atoms for the parameter"); + append_interaction(il, type, parm.atoms()); } } } @@ -573,7 +570,7 @@ void convertInteractionTypeParameters(int atnr, gmx::ArrayRef plist = intermolecular_interactions->plist; - if (plist[i].nr > 0) + if (!plist[i].interactionTypes.empty()) { flags = interaction_function[i].flags; /* For intermolecular interactions we (currently) diff --git a/src/gromacs/gmxpreprocess/gen_ad.cpp b/src/gromacs/gmxpreprocess/gen_ad.cpp index 9cb2eb7ffa..98d5031625 100644 --- a/src/gromacs/gmxpreprocess/gen_ad.cpp +++ b/src/gromacs/gmxpreprocess/gen_ad.cpp @@ -45,6 +45,7 @@ #include #include +#include #include "gromacs/fileio/confio.h" #include "gromacs/gmxpreprocess/gpp_nextnb.h" @@ -63,289 +64,185 @@ #include "resall.h" #define DIHEDRAL_WAS_SET_IN_RTP 0 -static bool was_dihedral_set_in_rtp(const t_param *dih) +static bool was_dihedral_set_in_rtp(const InteractionType &dih) { - return dih->c[MAXFORCEPARAM-1] == DIHEDRAL_WAS_SET_IN_RTP; + // This is a bad way to check this, but I don't know how to make this better now. + gmx::ArrayRef forceParam = dih.forceParam(); + return forceParam[MAXFORCEPARAM-1] == DIHEDRAL_WAS_SET_IN_RTP; } -typedef bool (*peq)(t_param *p1, t_param *p2); +typedef bool (*peq)(const InteractionType &p1, const InteractionType &p2); -static int acomp(const void *a1, const void *a2) +static bool acomp(const InteractionType &a1, const InteractionType &a2) { - const t_param *p1, *p2; int ac; - p1 = static_cast(a1); - p2 = static_cast(a2); - if ((ac = (p1->aj()-p2->aj())) != 0) + if ((ac = (a1.aj()-a2.aj())) != 0) { - return ac; + return ac < 0; } - else if ((ac = (p1->ai()-p2->ai())) != 0) + else if ((ac = (a1.ai()-a2.ai())) != 0) { - return ac; + return ac < 0; } else { - return (p1->ak()-p2->ak()); + return (a1.ak() < a2.ak()); } } -static int pcomp(const void *a1, const void *a2) +static bool pcomp(const InteractionType &a1, const InteractionType &a2) { - const t_param *p1, *p2; int pc; - p1 = static_cast(a1); - p2 = static_cast(a2); - if ((pc = (p1->ai()-p2->ai())) != 0) + if ((pc = (a1.ai()-a2.ai())) != 0) { - return pc; + return pc < 0; } else { - return (p1->aj()-p2->aj()); + return (a1.aj() < a2.aj()); } } -static int dcomp(const void *d1, const void *d2) +static bool dcomp(const InteractionType &d1, const InteractionType &d2) { - const t_param *p1, *p2; int dc; - p1 = static_cast(d1); - p2 = static_cast(d2); /* First sort by J & K (the two central) atoms */ - if ((dc = (p1->aj()-p2->aj())) != 0) + if ((dc = (d1.aj()-d2.aj())) != 0) { - return dc; + return dc < 0; } - else if ((dc = (p1->ak()-p2->ak())) != 0) + else if ((dc = (d1.ak()-d2.ak())) != 0) { - return dc; + return dc < 0; } /* Then make sure to put rtp dihedrals before generated ones */ - else if (was_dihedral_set_in_rtp(p1) && - !was_dihedral_set_in_rtp(p2)) + else if (was_dihedral_set_in_rtp(d1) && + !was_dihedral_set_in_rtp(d2)) { - return -1; + return true; } - else if (!was_dihedral_set_in_rtp(p1) && - was_dihedral_set_in_rtp(p2)) + else if (!was_dihedral_set_in_rtp(d1) && + was_dihedral_set_in_rtp(d2)) { - return 1; + return false; } /* Then sort by I and J (two outer) atoms */ - else if ((dc = (p1->ai()-p2->ai())) != 0) + else if ((dc = (d1.ai()-d2.ai())) != 0) { - return dc; + return dc < 0; } - else if ((dc = (p1->al()-p2->al())) != 0) + else if ((dc = (d1.al()-d2.al())) != 0) { - return dc; + return dc < 0; } else { // AMBER force fields with type 9 dihedrals can reach here, where we sort on // the contents of the string that names the macro for the parameters. - return strcmp(p1->s, p2->s); + return std::lexicographical_compare(d1.interactionTypeName().begin(), + d1.interactionTypeName().end(), + d2.interactionTypeName().begin(), + d2.interactionTypeName().end()); } } -static bool is_dihedral_on_same_bond(t_param *p1, t_param *p2) +static bool is_dihedral_on_same_bond(const InteractionType &p1, const InteractionType &p2) { - return ((p1->aj() == p2->aj()) && (p1->ak() == p2->ak())) || - ((p1->aj() == p2->ak()) && (p1->ak() == p2->aj())); + return ((p1.aj() == p2.aj()) && (p1.ak() == p2.ak())) || + ((p1.aj() == p2.ak()) && (p1.ak() == p2.aj())); } -static bool preq(t_param *p1, t_param *p2) +static bool preq(const InteractionType &p1, const InteractionType &p2) { - return (p1->ai() == p2->ai()) && (p1->aj() == p2->aj()); + return (p1.ai() == p2.ai()) && (p1.aj() == p2.aj()); } -static void rm2par(t_param p[], int *np, peq eq) +static void rm2par(std::vector *p, peq eq) { - int *index, nind; - int i, j; - - if ((*np) == 0) + if (p->empty()) { return; } - snew(index, *np); - nind = 0; - index[nind++] = 0; - for (i = 1; (i < (*np)); i++) - { - if (!eq(&p[i], &p[i-1])) - { - index[nind++] = i; - } - } - /* Index now holds pointers to all the non-equal params, - * this only works when p is sorted of course - */ - for (i = 0; (i < nind); i++) + for (auto param = p->begin() + 1; param != p->end(); ) { - for (j = 0; (j < MAXATOMLIST); j++) + auto prev = param - 1; + if (eq(*param, *prev)) { - p[i].a[j] = p[index[i]].a[j]; + param = p->erase(param); } - for (j = 0; (j < MAXFORCEPARAM); j++) - { - p[i].c[j] = p[index[i]].c[j]; - } - if (p[index[i]].a[0] == p[index[i]].a[1]) - { - strcpy(p[i].s, ""); - } - else if (index[i] > i) + else { - /* Copy the string only if it comes from somewhere else - * otherwise we will end up copying a random (newly freed) pointer. - * Since the index is sorted we only have to test for index[i] > i. - */ - strcpy(p[i].s, p[index[i]].s); + ++param; } } - (*np) = nind; - - sfree(index); } -static void cppar(t_param p[], int np, gmx::ArrayRef plist, int ftype) +static void cppar(gmx::ArrayRef types, + gmx::ArrayRef plist, + int ftype) { - int i, j, nral, nrfp; - - InteractionTypeParameters *ps = &plist[ftype]; - nral = NRAL(ftype); - nrfp = NRFP(ftype); - /* Keep old stuff */ - pr_alloc(np, ps); - for (i = 0; (i < np); i++) - { - for (j = 0; (j < nral); j++) - { - ps->param[ps->nr].a[j] = p[i].a[j]; - } - for (j = 0; (j < nrfp); j++) - { - ps->param[ps->nr].c[j] = p[i].c[j]; - } - for (j = 0; (j < MAXSLEN); j++) - { - ps->param[ps->nr].s[j] = p[i].s[j]; - } - ps->nr++; - } -} - -static void cpparam(t_param *dest, t_param *src) -{ - int j; - - for (j = 0; (j < MAXATOMLIST); j++) - { - dest->a[j] = src->a[j]; - } - for (j = 0; (j < MAXFORCEPARAM); j++) - { - dest->c[j] = src->c[j]; - } - for (j = 0; (j < MAXSLEN); j++) - { - dest->s[j] = src->s[j]; - } -} - -static void set_p(t_param *p, const int ai[4], const real *c, const char *s) -{ - int j; - - for (j = 0; (j < 4); j++) + for (const auto &type : types) { - p->a[j] = ai[j]; + plist[ftype].interactionTypes.push_back(type); } - for (j = 0; (j < MAXFORCEPARAM); j++) - { - if (c) - { - p->c[j] = c[j]; - } - else - { - p->c[j] = NOTSET; - } - } - - set_p_string(p, s); } -static int idcomp(const void *a, const void *b) +static bool idcomp(const InteractionType &a, const InteractionType &b) { - const t_param *pa, *pb; int d; - pa = static_cast(a); - pb = static_cast(b); - if ((d = (pa->a[0]-pb->a[0])) != 0) + if ((d = (a.ai()-b.ai())) != 0) { - return d; + return d < 0; } - else if ((d = (pa->a[3]-pb->a[3])) != 0) + else if ((d = (a.al()-b.al())) != 0) { - return d; + return d < 0; } - else if ((d = (pa->a[1]-pb->a[1])) != 0) + else if ((d = (a.aj()-b.aj())) != 0) { - return d; + return d < 0; } else { - return (pa->a[2]-pb->a[2]); + return (a.ak() < b.ak()); } } -static void sort_id(int nr, t_param ps[]) +static void sort_id(gmx::ArrayRef ps) { - int i, tmp; - - /* First swap order of atoms around if necessary */ - for (i = 0; (i < nr); i++) + if (ps.size() > 1) { - if (ps[i].a[3] < ps[i].a[0]) + for (auto &parm : ps) { - tmp = ps[i].a[3]; ps[i].a[3] = ps[i].a[0]; ps[i].a[0] = tmp; - tmp = ps[i].a[2]; ps[i].a[2] = ps[i].a[1]; ps[i].a[1] = tmp; + parm.sortAtomIds(); } - } - /* Now sort it */ - if (nr > 1) - { - qsort(ps, nr, static_cast(sizeof(ps[0])), idcomp); + std::sort(ps.begin(), ps.end(), idcomp); } } -static int n_hydro(const int a[], char ***atomname) +static int n_hydro(gmx::ArrayRef a, char ***atomname) { - int i, nh = 0; - char c0, c1, *aname; + int nh = 0; - for (i = 0; (i < 4); i += 3) + for (auto atom = a.begin(); atom < a.end(); atom += 3) { - aname = *atomname[a[i]]; - c0 = toupper(aname[0]); + const char *aname = *atomname[*atom]; + char c0 = toupper(aname[0]); if (c0 == 'H') { nh++; } else if ((static_cast(strlen(aname)) > 1) && (c0 >= '0') && (c0 <= '9')) { - c1 = toupper(aname[1]); + char c1 = toupper(aname[1]); if (c1 == 'H') { nh++; @@ -357,61 +254,56 @@ static int n_hydro(const int a[], char ***atomname) /* Clean up the dihedrals (both generated and read from the .rtp * file). */ -static void clean_dih(t_param *dih, int *ndih, t_param improper[], int nimproper, - t_atoms *atoms, bool bKeepAllGeneratedDihedrals, - bool bRemoveDihedralIfWithImproper) +static std::vector clean_dih(gmx::ArrayRef dih, + gmx::ArrayRef improper, + t_atoms *atoms, bool bKeepAllGeneratedDihedrals, + bool bRemoveDihedralIfWithImproper) { - int i, j, k, l; - int *index, nind; - /* Construct the list of the indices of the dihedrals * (i.e. generated or read) that might be kept. */ - snew(index, *ndih+1); + std::vector < std::pair < InteractionType, int>> newDihedrals; if (bKeepAllGeneratedDihedrals) { fprintf(stderr, "Keeping all generated dihedrals\n"); - nind = *ndih; - for (i = 0; i < nind; i++) + int i = 0; + for (const auto &dihedral : dih) { - index[i] = i; + newDihedrals.emplace_back(std::pair(dihedral, i++)); } - index[nind] = *ndih; } else { - nind = 0; /* Check if generated dihedral i should be removed. The * dihedrals have been sorted by dcomp() above, so all those * on the same two central atoms are together, with those from * the .rtp file preceding those that were automatically * generated. We remove the latter if the former exist. */ - for (i = 0; i < *ndih; i++) + int i = 0; + for (auto dihedral = dih.begin(); dihedral != dih.end(); dihedral++) { /* Keep the dihedrals that were defined in the .rtp file, * and the dihedrals that were generated and different * from the last one (whether it was generated or not). */ - if (was_dihedral_set_in_rtp(&dih[i]) || - 0 == i || - !is_dihedral_on_same_bond(&dih[i], &dih[i-1])) + if (was_dihedral_set_in_rtp(*dihedral) || + dihedral == dih.begin() || + !is_dihedral_on_same_bond(*dihedral, *(dihedral-1))) { - index[nind++] = i; + newDihedrals.emplace_back(std::pair(*dihedral, i++)); } } - index[nind] = *ndih; } - - k = 0; - for (i = 0; i < nind; i++) + int k = 0; + for (auto dihedral = newDihedrals.begin(); dihedral != newDihedrals.end(); ) { - bool bWasSetInRTP = was_dihedral_set_in_rtp(&dih[index[i]]); - bool bKeep = TRUE; + bool bWasSetInRTP = was_dihedral_set_in_rtp(dihedral->first); + bool bKeep = true; if (!bWasSetInRTP && bRemoveDihedralIfWithImproper) { /* Remove the dihedral if there is an improper on the same * bond. */ - for (j = 0; j < nimproper && bKeep; j++) + for (auto imp = improper.begin(); imp != improper.end() && bKeep; ++imp) { - bKeep = !is_dihedral_on_same_bond(&dih[index[i]], &improper[j]); + bKeep = !is_dihedral_on_same_bond(dihedral->first, *imp); } } @@ -424,17 +316,18 @@ static void clean_dih(t_param *dih, int *ndih, t_param improper[], int nimproper * index[]. However, their parameters are still present, * and l is looping over this dihedral and all of its * pruned siblings. */ - int bestl = index[i]; + int bestl = dihedral->second; if (!bKeepAllGeneratedDihedrals && !bWasSetInRTP) { /* Minimum number of hydrogens for i and l atoms */ int minh = 2; - for (l = index[i]; - (l < index[i+1] && - is_dihedral_on_same_bond(&dih[index[i]], &dih[l])); + int next = dihedral->second + 1; + for (int l = dihedral->second; + (l < next && + is_dihedral_on_same_bond(dihedral->first, dih[l])); l++) { - int nh = n_hydro(dih[l].a, atoms->atomname); + int nh = n_hydro(dih[l].atoms(), atoms->atomname); if (nh < minh) { minh = nh; @@ -445,37 +338,36 @@ static void clean_dih(t_param *dih, int *ndih, t_param improper[], int nimproper break; } } + dihedral->first = dih[bestl]; } - if (k != bestl) + if (k == bestl) { - cpparam(&dih[k], &dih[bestl]); + ++dihedral; } k++; } + else + { + dihedral = newDihedrals.erase(dihedral); + } } - - for (i = k; i < *ndih; i++) + std::vector finalDihedrals; + finalDihedrals.reserve(newDihedrals.size()); + for (const auto ¶m : newDihedrals) { - strcpy(dih[i].s, ""); + finalDihedrals.emplace_back(param.first); } - *ndih = k; - - sfree(index); + return finalDihedrals; } -static int get_impropers(t_atoms *atoms, gmx::ArrayRef globalPatches, t_param **improper, - bool bAllowMissing) +static std::vector get_impropers(t_atoms *atoms, + gmx::ArrayRef globalPatches, + bool bAllowMissing) { - int nimproper, start, ninc, nalloc; - int ai[MAXATOMLIST]; - - ninc = 500; - nalloc = ninc; - snew(*improper, nalloc); + std::vector improper; /* Add all the impropers from the residue database to the list. */ - nimproper = 0; - start = 0; + int start = 0; if (!globalPatches.empty()) { for (int i = 0; (i < atoms->nres); i++) @@ -483,27 +375,26 @@ static int get_impropers(t_atoms *atoms, gmx::ArrayRef gl BondedInteractionList *impropers = &globalPatches[i].rb[ebtsIDIHS]; for (const auto &bondeds : impropers->b) { - bool bStop = false; + bool bStop = false; + std::vector ai; for (int k = 0; (k < 4) && !bStop; k++) { - ai[k] = search_atom(bondeds.a[k].c_str(), start, - atoms, - "improper", bAllowMissing); - if (ai[k] == -1) + int entry = search_atom(bondeds.a[k].c_str(), start, + atoms, "improper", bAllowMissing); + + if (entry != -1) + { + ai.emplace_back(entry); + } + else { bStop = true; } } if (!bStop) { - if (nimproper == nalloc) - { - nalloc += ninc; - srenew(*improper, nalloc); - } /* Not broken out */ - set_p(&((*improper)[nimproper]), ai, nullptr, bondeds.s.c_str()); - nimproper++; + improper.emplace_back(InteractionType(ai, {}, bondeds.s)); } } while ((start < atoms->nr) && (atoms->atom[start].resind == i)) @@ -513,7 +404,7 @@ static int get_impropers(t_atoms *atoms, gmx::ArrayRef gl } } - return nimproper; + return improper; } static int nb_dist(t_nextnb *nnb, int ai, int aj) @@ -548,27 +439,25 @@ static bool is_hydro(t_atoms *atoms, int ai) return ((*(atoms->atomname[ai]))[0] == 'H'); } -static void get_atomnames_min(int n, char **anm, - int resind, t_atoms *atoms, const int *a) +static void get_atomnames_min(int n, gmx::ArrayRef anm, + int resind, t_atoms *atoms, gmx::ArrayRef a) { - int m; - /* Assume ascending residue numbering */ - for (m = 0; m < n; m++) + for (int m = 0; m < n; m++) { if (atoms->atom[a[m]].resind < resind) { - strcpy(anm[m], "-"); + anm[m] = "-"; } else if (atoms->atom[a[m]].resind > resind) { - strcpy(anm[m], "+"); + anm[m] = "+"; } else { - strcpy(anm[m], ""); + anm[m] = ""; } - strcat(anm[m], *(atoms->atomname[a[m]])); + anm[m].append(*(atoms->atomname[a[m]])); } } @@ -728,29 +617,16 @@ void gen_pad(t_nextnb *nnb, t_atoms *atoms, gmx::ArrayRef plist, t_excls excls[], gmx::ArrayRef globalPatches, bool bAllowMissing) { - t_param *ang, *dih, *pai, *improper; - char **anm; - int ninc, maxang, maxdih, maxpai; - int nang, ndih, npai, nimproper; - /* These are the angles, dihedrals and pairs that we generate * from the bonds. The ones that are already there from the rtp file * will be retained. */ - nang = 0; - npai = 0; - ndih = 0; - ninc = 500; - maxang = maxdih = maxpai = ninc; - snew(ang, maxang); - snew(dih, maxdih); - snew(pai, maxpai); - - snew(anm, 4); - for (int i = 0; i < 4; i++) - { - snew(anm[i], 12); - } + std::vector ang; + std::vector dih; + std::vector pai; + std::vector improper; + + std::array anm; if (!globalPatches.empty()) { @@ -786,31 +662,22 @@ void gen_pad(t_nextnb *nnb, t_atoms *atoms, gmx::ArrayRef atomNumbers = {i, j1, k1}; + std::string name; if (!globalPatches.empty()) { - int minres = atoms->atom[ang[nang].a[0]].resind; + int minres = atoms->atom[i].resind; int maxres = minres; for (int m = 1; m < 3; m++) { - minres = std::min(minres, atoms->atom[ang[nang].a[m]].resind); - maxres = std::max(maxres, atoms->atom[ang[nang].a[m]].resind); + minres = std::min(minres, atoms->atom[atomNumbers[m]].resind); + maxres = std::max(maxres, atoms->atom[atomNumbers[m]].resind); } int res = 2*minres-maxres; do { res += maxres-minres; - get_atomnames_min(3, anm, res, atoms, ang[nang].a); + get_atomnames_min(3, anm, res, atoms, atomNumbers); BondedInteractionList *hbang = &globalPatches[res].rb[ebtsANGLES]; for (auto &bondeds : hbang->b) { @@ -825,7 +692,7 @@ void gen_pad(t_nextnb *nnb, t_atoms *atoms, gmx::ArrayRefa[k1][1][l]; if ((l1 != i) && (l1 != j1)) { - if (ndih == maxdih) - { - maxdih += ninc; - srenew(dih, maxdih); - } - dih[ndih].ai() = i; - dih[ndih].aj() = j1; - dih[ndih].ak() = k1; - dih[ndih].al() = l1; - for (int m = 0; m < MAXFORCEPARAM; m++) - { - dih[ndih].c[m] = NOTSET; - } - set_p_string(&(dih[ndih]), ""); - int nFound = 0; + std::vector atomNumbers = {i, j1, k1, l1}; + std::string name; + int nFound = 0; if (!globalPatches.empty()) { - int minres = atoms->atom[dih[ndih].a[0]].resind; + int minres = atoms->atom[i].resind; int maxres = minres; for (int m = 1; m < 4; m++) { - minres = std::min(minres, atoms->atom[dih[ndih].a[m]].resind); - maxres = std::max(maxres, atoms->atom[dih[ndih].a[m]].resind); + minres = std::min( + minres, + atoms->atom[atomNumbers[m]].resind); + maxres = std::max( + maxres, + atoms->atom[atomNumbers[m]].resind); } int res = 2*minres-maxres; do { res += maxres-minres; - get_atomnames_min(4, anm, res, atoms, dih[ndih].a); + get_atomnames_min(4, anm, res, atoms, atomNumbers); BondedInteractionList *hbdih = &globalPatches[res].rb[ebtsPDIHS]; for (auto &bondeds : hbdih->b) { @@ -889,32 +748,16 @@ void gen_pad(t_nextnb *nnb, t_atoms *atoms, gmx::ArrayRef atoms = {i, j1, k1, l1}; + dih.push_back(InteractionType(atoms, {}, "")); } int nbd = nb_dist(nnb, i, l1); @@ -954,17 +784,8 @@ void gen_pad(t_nextnb *nnb, t_atoms *atoms, gmx::ArrayRef atoms = {i1, i2}; + pai.push_back(InteractionType(atoms, {}, "")); } } } @@ -995,12 +816,8 @@ void gen_pad(t_nextnb *nnb, t_atoms *atoms, gmx::ArrayRef atomNumbers; + bool bFound = true; for (int k = 0; k < 3 && bFound; k++) { const char *p = bondeds.a[k].c_str(); @@ -1015,18 +832,15 @@ void gen_pad(t_nextnb *nnb, t_atoms *atoms, gmx::ArrayRef atomNumbers; + bool bFound = true; for (int k = 0; k < 4 && bFound; k++) { const char *p = bondeds.a[k].c_str(); @@ -1060,83 +870,68 @@ void gen_pad(t_nextnb *nnb, t_atoms *atoms, gmx::ArrayRef 1) + if (ang.size() > 1) { - qsort(ang, nang, static_cast(sizeof(ang[0])), acomp); + std::sort(ang.begin(), ang.end(), acomp); } /* Sort dihedrals with respect to j-k-i-l (middle atoms first) */ - if (ndih > 1) + if (dih.size() > 1) { - qsort(dih, ndih, static_cast(sizeof(dih[0])), dcomp); + std::sort(dih.begin(), dih.end(), dcomp); } /* Sort the pairs */ - if (npai > 1) + if (pai.size() > 1) { - qsort(pai, npai, static_cast(sizeof(pai[0])), pcomp); + std::sort(pai.begin(), pai.end(), pcomp); } - if (npai > 0) + if (!pai.empty()) { /* Remove doubles, could occur in 6-rings, such as phenyls, maybe one does not want this when fudgeQQ < 1. */ - fprintf(stderr, "Before cleaning: %d pairs\n", npai); - rm2par(pai, &npai, preq); + fprintf(stderr, "Before cleaning: %zu pairs\n", pai.size()); + rm2par(&pai, preq); } /* Get the impropers from the database */ - nimproper = get_impropers(atoms, globalPatches, &improper, bAllowMissing); + improper = get_impropers(atoms, globalPatches, bAllowMissing); /* Sort the impropers */ - sort_id(nimproper, improper); + sort_id(improper); - if (ndih > 0) + if (!dih.empty()) { - fprintf(stderr, "Before cleaning: %d dihedrals\n", ndih); - clean_dih(dih, &ndih, improper, nimproper, atoms, - rtpFFDB[0].bKeepAllGeneratedDihedrals, - rtpFFDB[0].bRemoveDihedralIfWithImproper); + fprintf(stderr, "Before cleaning: %zu dihedrals\n", dih.size()); + dih = clean_dih(dih, improper, atoms, + rtpFFDB[0].bKeepAllGeneratedDihedrals, + rtpFFDB[0].bRemoveDihedralIfWithImproper); } /* Now we have unique lists of angles and dihedrals * Copy them into the destination struct */ - cppar(ang, nang, plist, F_ANGLES); - cppar(dih, ndih, plist, F_PDIHS); - cppar(improper, nimproper, plist, F_IDIHS); - cppar(pai, npai, plist, F_LJ14); + cppar(ang, plist, F_ANGLES); + cppar(dih, plist, F_PDIHS); + cppar(improper, plist, F_IDIHS); + cppar(pai, plist, F_LJ14); /* Remove all exclusions which are within nrexcl */ clean_excls(nnb, rtpFFDB[0].nrexcl, excls); - for (int i = 0; i < 4; i++) - { - sfree(anm[i]); - } - sfree(anm); - - sfree(ang); - sfree(dih); - sfree(improper); - sfree(pai); } diff --git a/src/gromacs/gmxpreprocess/gen_vsite.cpp b/src/gromacs/gmxpreprocess/gen_vsite.cpp index 1a50d6b100..b683cccd14 100644 --- a/src/gromacs/gmxpreprocess/gen_vsite.cpp +++ b/src/gromacs/gmxpreprocess/gen_vsite.cpp @@ -477,15 +477,15 @@ static void count_bonds(int atom, InteractionTypeParameters *psb, char ***atomna /* find heavy atom bound to this hydrogen */ heavy = NOTSET; - for (int i = 0; (i < psb->nr) && (heavy == NOTSET); i++) + for (auto parm = psb->interactionTypes.begin(); (parm != psb->interactionTypes.end()) && (heavy == NOTSET); parm++) { - if (psb->param[i].ai() == atom) + if (parm->ai() == atom) { - heavy = psb->param[i].aj(); + heavy = parm->aj(); } - else if (psb->param[i].aj() == atom) + else if (parm->aj() == atom) { - heavy = psb->param[i].ai(); + heavy = parm->ai(); } } if (heavy == NOTSET) @@ -497,15 +497,15 @@ static void count_bonds(int atom, InteractionTypeParameters *psb, char ***atomna nrb = 0; nrH = 0; nrhv = 0; - for (int i = 0; i < psb->nr; i++) + for (const auto &parm : psb->interactionTypes) { - if (psb->param[i].ai() == heavy) + if (parm.ai() == heavy) { - other = psb->param[i].aj(); + other = parm.aj(); } - else if (psb->param[i].aj() == heavy) + else if (parm.aj() == heavy) { - other = psb->param[i].ai(); + other = parm.ai(); } if (other != NOTSET) { @@ -669,15 +669,16 @@ static void add_vsites(gmx::ArrayRef plist, int vsite { /* find more heavy atoms */ other = moreheavy = NOTSET; - for (int j = 0; (j < plist[F_BONDS].nr) && (moreheavy == NOTSET); j++) + for (auto parm = plist[F_BONDS].interactionTypes.begin(); + (parm != plist[F_BONDS].interactionTypes.end()) && (moreheavy == NOTSET); parm++) { - if (plist[F_BONDS].param[j].ai() == heavies[0]) + if (parm->ai() == heavies[0]) { - other = plist[F_BONDS].param[j].aj(); + other = parm->aj(); } - else if (plist[F_BONDS].param[j].aj() == heavies[0]) + else if (parm->aj() == heavies[0]) { - other = plist[F_BONDS].param[j].ai(); + other = parm->ai(); } if ( (other != NOTSET) && (other != Heavy) ) { @@ -1568,7 +1569,6 @@ void do_vsites(gmx::ArrayRef rtpFFDB, PreprocessingAtom char tpname[32], nexttpname[32]; int *o2n, *newvsite_type, *newcgnr, ats[MAXATOMSPERRESIDUE]; t_atom *newatom; - InteractionTypeParameters *params; char ***newatomname; char *resnm = nullptr; int cmplength; @@ -2138,48 +2138,47 @@ void do_vsites(gmx::ArrayRef rtpFFDB, PreprocessingAtom InteractionTypeParameters *params = &(plist[ftype]); if (debug) { - fprintf(debug, "Renumbering %d %s\n", params->nr, + fprintf(debug, "Renumbering %zu %s\n", params->size(), interaction_function[ftype].longname); } - for (int i = 0; i < params->nr; i++) + /* Horrible hacks needed here to get this to work */ + for (auto parm = params->interactionTypes.begin(); parm != params->interactionTypes.end(); parm++) { + gmx::ArrayRef atomNumbers(parm->atoms()); + std::vector newAtomNumber; for (int j = 0; j < NRAL(ftype); j++) { - if (params->param[i].a[j] >= add_shift) + if (atomNumbers[j] >= add_shift) { if (debug) { - fprintf(debug, " [%d -> %d]", params->param[i].a[j], - params->param[i].a[j]-add_shift); + fprintf(debug, " [%d -> %d]", atomNumbers[j], + atomNumbers[j]-add_shift); } - params->param[i].a[j] = params->param[i].a[j]-add_shift; + newAtomNumber.emplace_back(atomNumbers[j]-add_shift); } else { if (debug) { - fprintf(debug, " [%d -> %d]", params->param[i].a[j], - o2n[params->param[i].a[j]]); + fprintf(debug, " [%d -> %d]", atomNumbers[j], + o2n[atomNumbers[j]]); } - params->param[i].a[j] = o2n[params->param[i].a[j]]; + newAtomNumber.emplace_back(o2n[atomNumbers[j]]); } } + *parm = InteractionType(newAtomNumber, parm->forceParam(), parm->interactionTypeName()); if (debug) { fprintf(debug, "\n"); } } } - /* now check if atoms in the added constraints are in increasing order */ - params = &(plist[F_CONSTRNC]); - for (int i = 0; i < params->nr; i++) + /* sort constraint parameters */ + InteractionTypeParameters *params = &(plist[F_CONSTRNC]); + for (auto &type : params->interactionTypes) { - if (params->param[i].ai() > params->param[i].aj()) - { - int j = params->param[i].aj(); - params->param[i].aj() = params->param[i].ai(); - params->param[i].ai() = j; - } + type.sortAtomIds(); } /* clean up */ @@ -2188,7 +2187,7 @@ void do_vsites(gmx::ArrayRef rtpFFDB, PreprocessingAtom /* tell the user what we did */ fprintf(stderr, "Marked %d virtual sites\n", nvsite); fprintf(stderr, "Added %d dummy masses\n", nadd); - fprintf(stderr, "Added %d new constraints\n", plist[F_CONSTRNC].nr); + fprintf(stderr, "Added %zu new constraints\n", plist[F_CONSTRNC].size()); } void do_h_mass(InteractionTypeParameters *psb, int vsite_type[], t_atoms *at, real mHmult, @@ -2202,18 +2201,18 @@ void do_h_mass(InteractionTypeParameters *psb, int vsite_type[], t_atoms *at, re { /* find bonded heavy atom */ int a = NOTSET; - for (int j = 0; (j < psb->nr) && (a == NOTSET); j++) + for (auto parm = psb->interactionTypes.begin(); (parm != psb->interactionTypes.end()) && (a == NOTSET); parm++) { /* if other atom is not a virtual site, it is the one we want */ - if ( (psb->param[j].ai() == i) && - !is_vsite(vsite_type[psb->param[j].aj()]) ) + if ( (parm->ai() == i) && + !is_vsite(vsite_type[parm->aj()]) ) { - a = psb->param[j].aj(); + a = parm->aj(); } - else if ( (psb->param[j].aj() == i) && - !is_vsite(vsite_type[psb->param[j].ai()]) ) + else if ( (parm->aj() == i) && + !is_vsite(vsite_type[parm->ai()]) ) { - a = psb->param[j].ai(); + a = parm->ai(); } } if (a == NOTSET) diff --git a/src/gromacs/gmxpreprocess/gpp_atomtype.cpp b/src/gromacs/gmxpreprocess/gpp_atomtype.cpp index 8a478dd0c5..df3781d39d 100644 --- a/src/gromacs/gmxpreprocess/gpp_atomtype.cpp +++ b/src/gromacs/gmxpreprocess/gpp_atomtype.cpp @@ -57,28 +57,27 @@ struct AtomTypeData { - //! 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), + //! Explicit constructor. + AtomTypeData(const t_atom &a, + char **name, + const InteractionType &nb, + const int bondAtomType, + const int atomNumber) : + atom_(a), name_(name), nb_(nb), 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_; + t_atom atom_; + //! Atom name. + char **name_; //! Nonbonded data. - t_param nb_; + InteractionType nb_; //! Bonded atomtype for the type. - int bondAtomType_; + int bondAtomType_; //! Atom number for the atom type. - int atomNumber_; + int atomNumber_; }; class PreprocessingAtomTypes::Impl @@ -162,11 +161,12 @@ real PreprocessingAtomTypes::atomNonBondedParamFromAtomType(int nt, int param) c { return NOTSET; } + gmx::ArrayRef forceParam = impl_->types[nt].nb_.forceParam(); if ((param < 0) || (param >= MAXFORCEPARAM)) { return NOTSET; } - return impl_->types[nt].nb_.c[param]; + return forceParam[param]; } PreprocessingAtomTypes::PreprocessingAtomTypes() @@ -186,12 +186,12 @@ PreprocessingAtomTypes &PreprocessingAtomTypes::operator=(PreprocessingAtomTypes PreprocessingAtomTypes::~PreprocessingAtomTypes() {} -int PreprocessingAtomTypes::addType(t_symtab *tab, - const t_atom &a, - const char *name, - const t_param *nb, - int bondAtomType, - int atomNumber) +int PreprocessingAtomTypes::addType(t_symtab *tab, + const t_atom &a, + const char *name, + const InteractionType &nb, + int bondAtomType, + int atomNumber) { auto found = std::find_if(impl_->types.begin(), impl_->types.end(), [&name](const AtomTypeData &data) @@ -212,13 +212,13 @@ int PreprocessingAtomTypes::addType(t_symtab *tab, } } -int PreprocessingAtomTypes::setType(int nt, - t_symtab *tab, - const t_atom &a, - const char *name, - const t_param *nb, - int bondAtomType, - int atomNumber) +int PreprocessingAtomTypes::setType(int nt, + t_symtab *tab, + const t_atom &a, + const char *name, + const InteractionType &nb, + int bondAtomType, + int atomNumber) { if (!isSet(nt)) { @@ -227,7 +227,7 @@ int PreprocessingAtomTypes::setType(int nt, impl_->types[nt].atom_ = a; impl_->types[nt].name_ = put_symtab(tab, name); - cp_param(&impl_->types[nt].nb_, nb); + impl_->types[nt].nb_ = nb; impl_->types[nt].bondAtomType_ = bondAtomType; impl_->types[nt].atomNumber_ = atomNumber; @@ -249,19 +249,18 @@ void PreprocessingAtomTypes::printTypes(FILE * out) fprintf (out, "\n"); } -static int search_atomtypes(const PreprocessingAtomTypes *ga, - int *n, - gmx::ArrayRef typelist, - int thistype, - t_param param[], - int ftype) +static int search_atomtypes(const PreprocessingAtomTypes *ga, + int *n, + gmx::ArrayRef typelist, + int thistype, + gmx::ArrayRef interactionTypes, + int ftype) { - int i, nn, nrfp, ntype; - - nn = *n; - nrfp = NRFP(ftype); - ntype = ga->size(); + int nn = *n; + int nrfp = NRFP(ftype); + int ntype = ga->size(); + int i; for (i = 0; (i < nn); i++) { if (typelist[i] == thistype) @@ -276,9 +275,11 @@ static int search_atomtypes(const PreprocessingAtomTypes *ga, for (int j = 0; j < ntype && bFound; j++) { /* Check nonbonded parameters */ - for (int k = 0; k < nrfp && bFound; k++) + gmx::ArrayRef forceParam1 = interactionTypes[ntype*typelist[i]+j].forceParam(); + gmx::ArrayRef forceParam2 = interactionTypes[ntype*thistype+j].forceParam(); + for (int k = 0; (k < nrfp) && bFound; k++) { - bFound = (param[ntype*typelist[i]+j].c[k] == param[ntype*thistype+j].c[k]); + bFound = forceParam1[k] == forceParam2[k]; } /* Check atomnumber */ @@ -330,7 +331,7 @@ void PreprocessingAtomTypes::renumberTypes(gmx::ArrayRef 0) + if (plist[F_LJ].size() > 0) { ftype = F_LJ; } @@ -351,10 +352,10 @@ void PreprocessingAtomTypes::renumberTypes(gmx::ArrayRefatom[i].type = search_atomtypes(this, &nat, typelist, atoms->atom[i].type, - plist[ftype].param, ftype); + plist[ftype].interactionTypes, ftype); atoms->atom[i].typeB = search_atomtypes(this, &nat, typelist, atoms->atom[i].typeB, - plist[ftype].param, ftype); + plist[ftype].interactionTypes, ftype); } } @@ -363,7 +364,7 @@ void PreprocessingAtomTypes::renumberTypes(gmx::ArrayRef= 0) { wall_atomtype[i] = search_atomtypes(this, &nat, typelist, wall_atomtype[i], - plist[ftype].param, ftype); + plist[ftype].interactionTypes, ftype); } } @@ -371,41 +372,24 @@ void PreprocessingAtomTypes::renumberTypes(gmx::ArrayRef nbsnew; - int nrfp = NRFP(ftype); - - int k = 0; for (int i = 0; (i < nat); i++) { int mi = typelist[i]; - for (int j = 0; (j < nat); j++, k++) + for (int j = 0; (j < nat); j++) { - int mj = typelist[j]; - for (int l = 0; (l < nrfp); l++) - { - nbsnew[k].c[l] = plist[ftype].param[ntype*mi+mj].c[l]; - } + int mj = typelist[j]; + const InteractionType &interactionType = plist[ftype].interactionTypes[ntype*mi+mj]; + nbsnew.emplace_back(interactionType.atoms(), interactionType.forceParam(), interactionType.interactionTypeName()); } new_types.push_back(impl_->types[mi]); } - int i; - for (i = 0; (i < nat*nat); i++) - { - for (int l = 0; (l < nrfp); l++) - { - plist[ftype].param[i].c[l] = nbsnew[i].c[l]; - } - } - plist[ftype].nr = i; mtop->ffparams.atnr = nat; - impl_->types = new_types; - - sfree(nbsnew); + impl_->types = new_types; + plist[ftype].interactionTypes = nbsnew; } void PreprocessingAtomTypes::copyTot_atomtypes(t_atomtypes *atomtypes) const diff --git a/src/gromacs/gmxpreprocess/gpp_atomtype.h b/src/gromacs/gmxpreprocess/gpp_atomtype.h index 57e0c95548..101653f38d 100644 --- a/src/gromacs/gmxpreprocess/gpp_atomtype.h +++ b/src/gromacs/gmxpreprocess/gpp_atomtype.h @@ -47,7 +47,7 @@ struct gmx_mtop_t; struct t_atom; struct t_atomtypes; -struct t_param; +class InteractionType; struct InteractionTypeParameters; struct t_symtab; @@ -181,13 +181,13 @@ class PreprocessingAtomTypes * \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); + int setType(int nt, + t_symtab *tab, + const t_atom &a, + const char *name, + const InteractionType &nb, + int bondAtomType, + int atomNumber); /*! \brief * Add new atom type to database. @@ -200,12 +200,12 @@ class PreprocessingAtomTypes * \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); + int addType(t_symtab *tab, + const t_atom &a, + const char *name, + const InteractionType &nb, + int bondAtomType, + int atomNumber); /*! \brief * Renumber existing atom type entries. diff --git a/src/gromacs/gmxpreprocess/gpp_nextnb.cpp b/src/gromacs/gmxpreprocess/gpp_nextnb.cpp index d793d94958..394c026e60 100644 --- a/src/gromacs/gmxpreprocess/gpp_nextnb.cpp +++ b/src/gromacs/gmxpreprocess/gpp_nextnb.cpp @@ -342,10 +342,11 @@ static void do_gen(int nrbonds, /* total number of bonds in s */ static void add_b(InteractionTypeParameters *bonds, int *nrf, sortable *s) { - for (int i = 0; (i < bonds->nr); i++) + int i = 0; + for (const auto &bond : bonds->interactionTypes) { - int ai = bonds->param[i].ai(); - int aj = bonds->param[i].aj(); + int ai = bond.ai(); + int aj = bond.aj(); if ((ai < 0) || (aj < 0)) { gmx_fatal(FARGS, "Impossible atom numbers in bond %d: ai=%d, aj=%d", @@ -356,6 +357,7 @@ static void add_b(InteractionTypeParameters *bonds, int *nrf, sortable *s) s[(*nrf)++].aj = aj; s[(*nrf)].aj = ai; s[(*nrf)++].ai = aj; + i++; } } @@ -370,7 +372,7 @@ void gen_nnb(t_nextnb *nnb, gmx::ArrayRef plist) if (IS_CHEMBOND(i)) { /* we need every bond twice (bidirectional) */ - nrbonds += 2*plist[i].nr; + nrbonds += 2*plist[i].size(); } } diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index 1d6430be8b..8b01bd5f25 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -105,9 +105,129 @@ #include "gromacs/utility/smalloc.h" #include "gromacs/utility/snprintf.h" +/* TODO The implementation details should move to their own source file. */ +InteractionType::InteractionType(gmx::ArrayRef atoms, + gmx::ArrayRef params, + const std::string &name) + : atoms_(atoms.begin(), atoms.end()), + interactionTypeName_(name) +{ + GMX_RELEASE_ASSERT(params.size() <= forceParam_.size(), + gmx::formatString("Cannot have more parameters than the maximum number possible (%d)", MAXFORCEPARAM).c_str()); + auto forceParamIt = forceParam_.begin(); + for (const auto param : params) + { + *forceParamIt++ = param; + } + while (forceParamIt != forceParam_.end()) + { + *forceParamIt++ = NOTSET; + } +} + +const int &InteractionType::ai() const +{ + GMX_RELEASE_ASSERT(!atoms_.empty(), "Need to have at least one atom set"); + return atoms_[0]; +} + +const int &InteractionType::aj() const +{ + GMX_RELEASE_ASSERT(atoms_.size() > 1, "Need to have at least two atoms set"); + return atoms_[1]; +} + +const int &InteractionType::ak() const +{ + GMX_RELEASE_ASSERT(atoms_.size() > 2, "Need to have at least three atoms set"); + return atoms_[2]; +} + +const int &InteractionType::al() const +{ + GMX_RELEASE_ASSERT(atoms_.size() > 3, "Need to have at least four atoms set"); + return atoms_[3]; +} + +const int &InteractionType::am() const +{ + GMX_RELEASE_ASSERT(atoms_.size() > 4, "Need to have at least five atoms set"); + return atoms_[4]; +} + +const real &InteractionType::c0() const +{ + return forceParam_[0]; +} + +const real &InteractionType::c1() const +{ + return forceParam_[1]; +} + +const real &InteractionType::c2() const +{ + return forceParam_[2]; +} + +const std::string &InteractionType::interactionTypeName() const +{ + return interactionTypeName_; +} + +void InteractionType::sortBondAtomIds() +{ + if (aj() < ai()) + { + std::swap(atoms_[0], atoms_[1]); + } +} + +void InteractionType::sortAngleAtomIds() +{ + if (ak() < ai()) + { + std::swap(atoms_[0], atoms_[2]); + } +} + +void InteractionType::sortDihedralAtomIds() +{ + if (al() < ai()) + { + std::swap(atoms_[0], atoms_[3]); + std::swap(atoms_[1], atoms_[2]); + } +} + +void InteractionType::sortAtomIds() +{ + if (isBond()) + { + sortBondAtomIds(); + } + else if (isAngle()) + { + sortAngleAtomIds(); + } + else if (isDihedral()) + { + sortDihedralAtomIds(); + } + else + { + GMX_THROW(gmx::InternalError("Can not sort parameters other than bonds, angles or dihedrals")); + } +}; + +void InteractionType::setForceParameter(int pos, real value) +{ + GMX_RELEASE_ASSERT(pos < MAXFORCEPARAM, "Can't set parameter beyond the maximum number of parameters"); + forceParam_[pos] = value; +} + void MoleculeInformation::initMolInfo() { - init_plist(plist); init_block(&cgs); init_block(&mols); init_blocka(&excls); @@ -117,7 +237,6 @@ void MoleculeInformation::initMolInfo() void MoleculeInformation::partialCleanUp() { done_block(&mols); - done_plist(plist); } void MoleculeInformation::fullCleanUp() @@ -125,7 +244,6 @@ void MoleculeInformation::fullCleanUp() done_atom (&atoms); done_block(&cgs); done_block(&mols); - done_plist(plist); } static int rm_interactions(int ifunc, gmx::ArrayRef mols) @@ -134,8 +252,8 @@ static int rm_interactions(int ifunc, gmx::ArrayRef mols) /* For all the molecule types */ for (auto &mol : mols) { - n += mol.plist[ifunc].nr; - mol.plist[ifunc].nr = 0; + n += mol.plist[ifunc].size(); + mol.plist[ifunc].interactionTypes.clear(); } return n; } @@ -433,7 +551,7 @@ static int nint_ftype(gmx_mtop_t *mtop, gmx::ArrayRef int nint = 0; for (const gmx_molblock_t &molb : mtop->molblock) { - nint += molb.nmol*mi[molb.type].plist[ftype].nr; + nint += molb.nmol*mi[molb.type].plist[ftype].size(); } return nint; @@ -836,7 +954,7 @@ static void read_posres(gmx_mtop_t *mtop, matrix box, invbox; int natoms, npbcdim = 0; char warn_buf[STRLEN]; - int a, i, ai, j, k, nat_molb; + int a, nat_molb; t_atom *atom; snew(top, 1); @@ -855,7 +973,7 @@ static void read_posres(gmx_mtop_t *mtop, if (rc_scaling != erscNO) { copy_mat(box, invbox); - for (j = npbcdim; j < DIM; j++) + for (int j = npbcdim; j < DIM; j++) { clear_rvec(invbox[j]); invbox[j][j] = 1; @@ -873,12 +991,12 @@ static void read_posres(gmx_mtop_t *mtop, nat_molb = molb.nmol*mtop->moltype[molb.type].atoms.nr; const InteractionTypeParameters *pr = &(molinfo[molb.type].plist[F_POSRES]); const InteractionTypeParameters *prfb = &(molinfo[molb.type].plist[F_FBPOSRES]); - if (pr->nr > 0 || prfb->nr > 0) + if (pr->size() > 0 || prfb->size() > 0) { atom = mtop->moltype[molb.type].atoms.atom; - for (i = 0; (i < pr->nr); i++) + for (const auto &restraint : pr->interactionTypes) { - ai = pr->param[i].ai(); + int ai = restraint.ai(); if (ai >= natoms) { gmx_fatal(FARGS, "Position restraint atom index (%d) in moltype '%s' is larger than number of atoms in %s (%d).\n", @@ -888,7 +1006,7 @@ static void read_posres(gmx_mtop_t *mtop, if (rc_scaling == erscCOM) { /* Determine the center of mass of the posres reference coordinates */ - for (j = 0; j < npbcdim; j++) + for (int j = 0; j < npbcdim; j++) { sum[j] += atom[ai].m*x[a+ai][j]; } @@ -896,9 +1014,9 @@ static void read_posres(gmx_mtop_t *mtop, } } /* Same for flat-bottomed posres, but do not count an atom twice for COM */ - for (i = 0; (i < prfb->nr); i++) + for (const auto &restraint : prfb->interactionTypes) { - ai = prfb->param[i].ai(); + int ai = restraint.ai(); if (ai >= natoms) { gmx_fatal(FARGS, "Position restraint atom index (%d) in moltype '%s' is larger than number of atoms in %s (%d).\n", @@ -907,7 +1025,7 @@ static void read_posres(gmx_mtop_t *mtop, if (rc_scaling == erscCOM && !hadAtom[ai]) { /* Determine the center of mass of the posres reference coordinates */ - for (j = 0; j < npbcdim; j++) + for (int j = 0; j < npbcdim; j++) { sum[j] += atom[ai].m*x[a+ai][j]; } @@ -917,7 +1035,7 @@ static void read_posres(gmx_mtop_t *mtop, if (!bTopB) { molb.posres_xA.resize(nat_molb); - for (i = 0; i < nat_molb; i++) + for (int i = 0; i < nat_molb; i++) { copy_rvec(x[a+i], molb.posres_xA[i]); } @@ -925,7 +1043,7 @@ static void read_posres(gmx_mtop_t *mtop, else { molb.posres_xB.resize(nat_molb); - for (i = 0; i < nat_molb; i++) + for (int i = 0; i < nat_molb; i++) { copy_rvec(x[a+i], molb.posres_xB[i]); } @@ -939,7 +1057,7 @@ static void read_posres(gmx_mtop_t *mtop, { gmx_fatal(FARGS, "The total mass of the position restraint atoms is 0"); } - for (j = 0; j < npbcdim; j++) + for (int j = 0; j < npbcdim; j++) { com[j] = sum[j]/totmass; } @@ -956,15 +1074,15 @@ static void read_posres(gmx_mtop_t *mtop, if (!molb.posres_xA.empty() || !molb.posres_xB.empty()) { std::vector &xp = (!bTopB ? molb.posres_xA : molb.posres_xB); - for (i = 0; i < nat_molb; i++) + for (int i = 0; i < nat_molb; i++) { - for (j = 0; j < npbcdim; j++) + for (int j = 0; j < npbcdim; j++) { if (rc_scaling == erscALL) { /* Convert from Cartesian to crystal coordinates */ xp[i][j] *= invbox[j][j]; - for (k = j+1; k < npbcdim; k++) + for (int k = j+1; k < npbcdim; k++) { xp[i][j] += invbox[k][j]*xp[i][k]; } @@ -982,10 +1100,10 @@ static void read_posres(gmx_mtop_t *mtop, if (rc_scaling == erscCOM) { /* Convert the COM from Cartesian to crystal coordinates */ - for (j = 0; j < npbcdim; j++) + for (int j = 0; j < npbcdim; j++) { com[j] *= invbox[j][j]; - for (k = j+1; k < npbcdim; k++) + for (int k = j+1; k < npbcdim; k++) { com[j] += invbox[k][j]*com[k]; } @@ -1207,7 +1325,7 @@ static int count_constraints(const gmx_mtop_t *mtop, gmx::ArrayRef mi, warninp *wi) { - int count, count_mol, i; + int count, count_mol; char buf[STRLEN]; count = 0; @@ -1216,15 +1334,15 @@ static int count_constraints(const gmx_mtop_t *mtop, count_mol = 0; gmx::ArrayRef plist = mi[molb.type].plist; - for (i = 0; i < F_NRE; i++) + for (int i = 0; i < F_NRE; i++) { if (i == F_SETTLE) { - count_mol += 3*plist[i].nr; + count_mol += 3*plist[i].size(); } else if (interaction_function[i].flags & IF_CONSTRAINT) { - count_mol += plist[i].nr; + count_mol += plist[i].size(); } } @@ -1797,7 +1915,6 @@ int gmx_grompp(int argc, char *argv[]) } std::array plist; - init_plist(plist); gmx_mtop_t sys; PreprocessingAtomTypes atypes; if (debug) @@ -2319,7 +2436,6 @@ int gmx_grompp(int argc, char *argv[]) sfree(opts->define); sfree(opts->include); sfree(opts); - done_plist(plist); for (auto &mol : mi) { // Some of the contents of molinfo have been stolen, so diff --git a/src/gromacs/gmxpreprocess/grompp_impl.h b/src/gromacs/gmxpreprocess/grompp_impl.h index e2b3c55fe6..eae3711ae5 100644 --- a/src/gromacs/gmxpreprocess/grompp_impl.h +++ b/src/gromacs/gmxpreprocess/grompp_impl.h @@ -40,71 +40,115 @@ #include +#include "gromacs/gmxpreprocess/notset.h" #include "gromacs/topology/atoms.h" #include "gromacs/topology/block.h" #include "gromacs/topology/idef.h" +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/basedefinitions.h" +#include "gromacs/utility/exceptions.h" #include "gromacs/utility/real.h" -#define MAXSLEN 32 - -struct t_param +/*! \libinternal \brief + * Describes a single parameter in a force field. + */ +class InteractionType { - int a[MAXATOMLIST]; /* The atom list (eg. bonds: particle */ - /* i = a[0] (ai), j = a[1] (aj)) */ - real c[MAXFORCEPARAM]; /* Force parameters (eg. b0 = c[0]) */ - char s[MAXSLEN]; /* A string (instead of parameters), * - * read from the .rtp file in pdb2gmx */ - const int &ai() const { return a[0]; } - int &ai() { return a[0]; } - const int &aj() const { return a[1]; } - int &aj() { return a[1]; } - const int &ak() const { return a[2]; } - int &ak() { return a[2]; } - const int &al() const { return a[3]; } - int &al() { return a[3]; } - const int &am() const { return a[4]; } - int &am() { return a[4]; } - - const real &c0() const { return c[0]; } - real &c0() { return c[0]; } - const real &c1() const { return c[1]; } - real &c1() { return c[1]; } - const real &c2() const { return c[2]; } - real &c2() { return c[2]; } + public: + //! Constructor that initializes vectors. + InteractionType(gmx::ArrayRef atoms, + gmx::ArrayRef params, + const std::string &name = ""); + /*!@{*/ + //! Access the individual elements set for the parameter. + const int &ai() const; + const int &aj() const; + const int &ak() const; + const int &al() const; + const int &am() const; + + const real &c0() const; + const real &c1() const; + const real &c2() const; + + const std::string &interactionTypeName() const; + /*!@}*/ + + /*! \brief Renumbers atom Ids. + * + * Enforces that ai() is less than the opposite terminal atom index, + * with the number depending on the interaction type. + */ + void sortAtomIds(); + + //! Set single force field parameter. + void setForceParameter(int pos, real value); + + //! View on all atoms numbers that are actually set. + gmx::ArrayRef atoms() { return atoms_; } + //! Const view on all atoms numbers that are actually set. + gmx::ArrayRef atoms() const { return atoms_; } + //! View on all of the force field parameters + gmx::ArrayRef forceParam() const { return forceParam_; } + //! View on all of the force field parameters + gmx::ArrayRef forceParam() { return forceParam_; } + + private: + //! Return if we have a bond parameter, means two atoms right now. + bool isBond() const { return atoms_.size() == 2; } + //! Return if we have an angle parameter, means three atoms right now. + bool isAngle() const { return atoms_.size() == 3; } + //! Return if we have a dihedral parameter, means four atoms right now. + bool isDihedral() const { return atoms_.size() == 4; } + //! Return if we have a cmap parameter, means five atoms right now. + bool isCmap() const { return atoms_.size() == 5; } + //! Enforce that atom id ai() is less than aj(). + void sortBondAtomIds(); + //! Enforce that atom id ai() is less than ak(). Does not change aj(). + void sortAngleAtomIds(); + /*! \brief Enforce order of atoms in dihedral. + * + * Changes atom order if needed to enforce that ai() is less than al(). + * If ai() and al() are swapped, aj() and ak() are swapped as well, + * independent of their previous order. + */ + void sortDihedralAtomIds(); + //! The atom list (eg. bonds: particle, i = atoms[0] (ai), j = atoms[1] (aj)) + std::vector atoms_; + //! Force parameters (eg. b0 = forceParam[0]) + std::array forceParam_; + //! Used with forcefields whose .rtp files name the interaction types (e.g. GROMOS), rather than look them up from the atom names. + std::string interactionTypeName_; }; /*! \libinternal \brief - * The parameters for a set of interactions of a given - * (unspecified) type found in the enumeration in ifunc.h. + * A set of interactions of a given type + * (found in the enumeration in ifunc.h), complete with + * atom indices and force field function parameters. * * This is used for containing the data obtained from the * lists of interactions of a given type in a [moleculetype] * topology file definition. - * - * \todo Remove manual memory management. */ struct InteractionTypeParameters { // NOLINT (clang-analyzer-optin.performance.Padding) - //! Number of parameters. - int nr = 0; - //! Maximum number of parameters. - int maxnr = 0; - //! The parameters of the interactions - t_param *param; + //! The different parameters in the system. + std::vector interactionTypes; //! CMAP grid spacing. - int cmakeGridSpacing = -1; + int cmakeGridSpacing = -1; //! Number of cmap angles. - int cmapAngles = -1; + int cmapAngles = -1; //! CMAP grid data. - std::vector cmap; + std::vector cmap; //! The five atomtypes followed by a number that identifies the type. - std::vector cmapAtomTypes; + std::vector cmapAtomTypes; + //! Number of parameters. + size_t size() const { return interactionTypes.size(); } //! Elements in cmap grid data. - int ncmap() const { return cmap.size(); } + int ncmap() const { return cmap.size(); } //! Number of elements in cmapAtomTypes. - int nct() const { return cmapAtomTypes.size(); } + int nct() const { return cmapAtomTypes.size(); } }; struct t_excls diff --git a/src/gromacs/gmxpreprocess/nm2type.cpp b/src/gromacs/gmxpreprocess/nm2type.cpp index 27c0a8a1c2..9a7513d9f0 100644 --- a/src/gromacs/gmxpreprocess/nm2type.cpp +++ b/src/gromacs/gmxpreprocess/nm2type.cpp @@ -194,14 +194,12 @@ static int match_str(const char *atom, const char *template_string) int nm2type(int nnm, t_nm2type nm2t[], t_symtab *tab, t_atoms *atoms, PreprocessingAtomTypes *atype, int *nbonds, InteractionTypeParameters *bonds) { - int cur = 0; + int cur = 0; #define prev (1-cur) - int nresolved, nb, maxbond, nqual[2][ematchNR]; - int *bbb, *n_mask, *m_mask, **match; - char *aname_i, *aname_n; - t_param *param; + int nresolved, nb, maxbond, nqual[2][ematchNR]; + int *bbb, *n_mask, *m_mask, **match; + char *aname_i, *aname_n; - snew(param, 1); maxbond = 0; for (int i = 0; (i < atoms->nr); i++) { @@ -225,10 +223,10 @@ int nm2type(int nnm, t_nm2type nm2t[], t_symtab *tab, t_atoms *atoms, { aname_i = *atoms->atomname[i]; nb = 0; - for (int j = 0; (j < bonds->nr); j++) + for (const auto &bond : bonds->interactionTypes) { - int ai = bonds->param[j].ai(); - int aj = bonds->param[j].aj(); + int ai = bond.ai(); + int aj = bond.aj(); if (ai == i) { bbb[nb++] = aj; @@ -338,7 +336,7 @@ int nm2type(int nnm, t_nm2type nm2t[], t_symtab *tab, t_atoms *atoms, { atoms->atom[i].qB = alpha; atoms->atom[i].m = atoms->atom[i].mB = mm; - k = atype->addType(tab, atoms->atom[i], type, param, + k = atype->addType(tab, atoms->atom[i], type, InteractionType({}, {}), atoms->atom[i].type, atomnr); } atoms->atom[i].type = k; @@ -358,7 +356,6 @@ int nm2type(int nnm, t_nm2type nm2t[], t_symtab *tab, t_atoms *atoms, sfree(bbb); sfree(n_mask); sfree(m_mask); - sfree(param); return nresolved; } diff --git a/src/gromacs/gmxpreprocess/pdb2top.cpp b/src/gromacs/gmxpreprocess/pdb2top.cpp index f2f24f03d4..bf7ea227ca 100644 --- a/src/gromacs/gmxpreprocess/pdb2top.cpp +++ b/src/gromacs/gmxpreprocess/pdb2top.cpp @@ -727,7 +727,7 @@ static void do_ssbonds(InteractionTypeParameters *ps, t_atoms *atoms, gmx_fatal(FARGS, "Trying to make impossible special bond (%s-%s)!", bond.firstAtom.c_str(), bond.secondAtom.c_str()); } - add_param(ps, ai, aj, nullptr, nullptr); + add_param(ps, ai, aj, {}, nullptr); } } @@ -780,7 +780,7 @@ static void at2bonds(InteractionTypeParameters *psb, gmx::ArrayRef(a); - pb = static_cast(b); + int d; - d = pa->a[0] - pb->a[0]; - if (d == 0) + if ((d = a.ai() - b.ai()) != 0) { - d = pa->a[1] - pb->a[1]; + return d < 0; } - if (d == 0) + else if ((d = a.aj() - b.aj()) != 0) { - return strlen(pb->s) - strlen(pa->s); + return d < 0; } else { - return d; + return a.interactionTypeName().length() > b.interactionTypeName().length(); } } static void clean_bonds(InteractionTypeParameters *ps) { - int i, j; - int a; - - if (ps->nr > 0) + if (ps->size() > 0) { - /* swap atomnumbers in bond if first larger than second: */ - for (i = 0; (i < ps->nr); i++) + /* Sort bonds */ + for (auto &bond : ps->interactionTypes) { - if (ps->param[i].a[1] < ps->param[i].a[0]) - { - a = ps->param[i].a[0]; - ps->param[i].a[0] = ps->param[i].a[1]; - ps->param[i].a[1] = a; - } + bond.sortAtomIds(); } - - /* Sort bonds */ - qsort(ps->param, ps->nr, static_cast(sizeof(ps->param[0])), pcompar); + std::sort(ps->interactionTypes.begin(), ps->interactionTypes.end(), pcompar); /* remove doubles, keep the first one always. */ - j = 1; - for (i = 1; (i < ps->nr); i++) + int oldNumber = ps->size(); + for (auto parm = ps->interactionTypes.begin() + 1; parm != ps->interactionTypes.end(); ) { - if ((ps->param[i].a[0] != ps->param[j-1].a[0]) || - (ps->param[i].a[1] != ps->param[j-1].a[1]) ) + auto prev = parm - 1; + if (parm->ai() == prev->ai() && + parm->aj() == prev->aj()) { - if (j != i) - { - cp_param(&(ps->param[j]), &(ps->param[i])); - } - j++; + parm = ps->interactionTypes.erase(parm); + } + else + { + ++parm; } } - fprintf(stderr, "Number of bonds was %d, now %d\n", ps->nr, j); - ps->nr = j; + fprintf(stderr, "Number of bonds was %d, now %zu\n", oldNumber, ps->size()); } else { @@ -1486,7 +1472,6 @@ void pdb2top(FILE *top_file, const char *posre_fn, const char *molname, int i, nmissat; int bts[ebtsNR]; - init_plist(plist); ResidueType rt; /* Make bonds */ @@ -1548,10 +1533,10 @@ void pdb2top(FILE *top_file, const char *posre_fn, const char *molname, if (bCmap) { gen_cmap(&(plist[F_CMAP]), usedPpResidues, atoms); - if (plist[F_CMAP].nr > 0) + if (plist[F_CMAP].size() > 0) { - fprintf(stderr, "There are %4d cmap torsion pairs\n", - plist[F_CMAP].nr); + fprintf(stderr, "There are %4zu cmap torsion pairs\n", + plist[F_CMAP].size()); } } @@ -1566,17 +1551,17 @@ void pdb2top(FILE *top_file, const char *posre_fn, const char *molname, /* clean_bonds(&(plist[F_BONDS]));*/ fprintf(stderr, - "There are %4d dihedrals, %4d impropers, %4d angles\n" - " %4d pairs, %4d bonds and %4d virtual sites\n", - plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr, - plist[F_LJ14].nr, plist[F_BONDS].nr, - plist[F_VSITE2].nr + - plist[F_VSITE3].nr + - plist[F_VSITE3FD].nr + - plist[F_VSITE3FAD].nr + - plist[F_VSITE3OUT].nr + - plist[F_VSITE4FD].nr + - plist[F_VSITE4FDN].nr ); + "There are %4zu dihedrals, %4zu impropers, %4zu angles\n" + " %4zu pairs, %4zu bonds and %4zu virtual sites\n", + plist[F_PDIHS].size(), plist[F_IDIHS].size(), plist[F_ANGLES].size(), + plist[F_LJ14].size(), plist[F_BONDS].size(), + plist[F_VSITE2].size() + + plist[F_VSITE3].size() + + plist[F_VSITE3FD].size() + + plist[F_VSITE3FAD].size() + + plist[F_VSITE3OUT].size() + + plist[F_VSITE4FD].size() + + plist[F_VSITE4FDN].size() ); print_sums(atoms, FALSE); @@ -1612,10 +1597,6 @@ void pdb2top(FILE *top_file, const char *posre_fn, const char *molname, /* we should clean up hb and restp here, but that is a *L*O*T* of work! */ sfree(cgnr); - for (i = 0; i < F_NRE; i++) - { - sfree(plist[i].param); - } for (i = 0; i < atoms->nr; i++) { sfree(excls[i].e); diff --git a/src/gromacs/gmxpreprocess/resall.cpp b/src/gromacs/gmxpreprocess/resall.cpp index c276afa5da..406a2db1c6 100644 --- a/src/gromacs/gmxpreprocess/resall.cpp +++ b/src/gromacs/gmxpreprocess/resall.cpp @@ -68,12 +68,10 @@ PreprocessingAtomTypes read_atype(const char *ffdir, t_symtab *tab) double m; int nratt = 0; t_atom *a; - t_param *nb; std::vector files = fflib_search_file_end(ffdir, ".atp", TRUE); snew(a, 1); - snew(nb, 1); - PreprocessingAtomTypes at; + PreprocessingAtomTypes at; for (const auto &filename : files) { @@ -94,7 +92,7 @@ PreprocessingAtomTypes read_atype(const char *ffdir, t_symtab *tab) if (sscanf(buf, "%s%lf", name, &m) == 2) { a->m = m; - at.addType(tab, *a, name, nb, 0, 0); + at.addType(tab, *a, name, InteractionType({}, {}), 0, 0); fprintf(stderr, "\rAtomtype %d", ++nratt); fflush(stderr); } diff --git a/src/gromacs/gmxpreprocess/tomorse.cpp b/src/gromacs/gmxpreprocess/tomorse.cpp index 934a773b0c..5c854a88a6 100644 --- a/src/gromacs/gmxpreprocess/tomorse.cpp +++ b/src/gromacs/gmxpreprocess/tomorse.cpp @@ -192,8 +192,6 @@ void convert_harmonics(gmx::ArrayRef mols, PreprocessingAto int n2m; t_2morse *t2m; - bool *bRemoveHarm; - /* First get the data */ t2m = read_dissociation_energies(&n2m); if (n2m <= 0) @@ -209,65 +207,42 @@ void convert_harmonics(gmx::ArrayRef mols, PreprocessingAto /* Check how many morse and harmonic BONDSs there are, increase size of * morse with the number of harmonics */ - int nrmorse = mol.plist[F_MORSE].nr; - for (int bb = 0; (bb < F_NRE); bb++) { if ((interaction_function[bb].flags & IF_BTYPE) && (bb != F_MORSE)) { - int nrharm = mol.plist[bb].nr; - pr_alloc(nrharm, &(mol.plist[F_MORSE])); - snew(bRemoveHarm, nrharm); + int nrharm = mol.plist[bb].size(); /* Now loop over the harmonics, trying to convert them */ - for (int j = 0; (j < nrharm); j++) + for (auto harmonic = mol.plist[bb].interactionTypes.begin(); + harmonic != mol.plist[bb].interactionTypes.end(); ) { - int ni = mol.plist[bb].param[j].ai(); - int nj = mol.plist[bb].param[j].aj(); + int ni = harmonic->ai(); + int nj = harmonic->aj(); real edis = search_e_diss(n2m, t2m, atype->atomNameFromAtomType(mol.atoms.atom[ni].type), atype->atomNameFromAtomType(mol.atoms.atom[nj].type)); if (edis != 0) { - bRemoveHarm[j] = true; - real b0 = mol.plist[bb].param[j].c[0]; - real kb = mol.plist[bb].param[j].c[1]; - real beta = std::sqrt(kb/(2*edis)); - mol.plist[F_MORSE].param[nrmorse].a[0] = ni; - mol.plist[F_MORSE].param[nrmorse].a[1] = nj; - mol.plist[F_MORSE].param[nrmorse].c[0] = b0; - mol.plist[F_MORSE].param[nrmorse].c[1] = edis; - mol.plist[F_MORSE].param[nrmorse].c[2] = beta; - nrmorse++; + real b0 = harmonic->c0(); + real kb = harmonic->c1(); + real beta = std::sqrt(kb/(2*edis)); + std::vector atoms = {ni, nj}; + std::vector forceParam = {b0, edis, beta}; + mol.plist[F_MORSE].interactionTypes.emplace_back( + InteractionType(atoms, forceParam)); + harmonic = mol.plist[bb].interactionTypes.erase(harmonic); } - } - mol.plist[F_MORSE].nr = nrmorse; - - /* Now remove the harmonics */ - int last = 0; - for (int j = 0; (j < nrharm); j++) - { - if (!bRemoveHarm[j]) + else { - /* Copy it to the last position */ - for (int k = 0; (k < MAXATOMLIST); k++) - { - mol.plist[bb].param[last].a[k] = - mol.plist[bb].param[j].a[k]; - } - for (int k = 0; (k < MAXFORCEPARAM); k++) - { - mol.plist[bb].param[last].c[k] = - mol.plist[bb].param[j].c[k]; - } - last++; + ++harmonic; } } - sfree(bRemoveHarm); + + int newHarmonics = mol.plist[bb].size(); fprintf(stderr, "Converted %d out of %d %s to morse bonds for mol %d\n", - nrharm-last, nrharm, interaction_function[bb].name, i); - mol.plist[bb].nr = last; + nrharm-newHarmonics, nrharm, interaction_function[bb].name, i); } } i++; diff --git a/src/gromacs/gmxpreprocess/topio.cpp b/src/gromacs/gmxpreprocess/topio.cpp index eb0beb201f..ca630939ff 100644 --- a/src/gromacs/gmxpreprocess/topio.cpp +++ b/src/gromacs/gmxpreprocess/topio.cpp @@ -88,13 +88,12 @@ static void gen_pairs(const InteractionTypeParameters &nbs, InteractionTypeParameters *pairs, real fudge, int comb) { real scaling; - int ntp = nbs.nr; + int ntp = nbs.size(); 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"); 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)) { @@ -102,17 +101,18 @@ static void gen_pairs(const InteractionTypeParameters &nbs, InteractionTypeParam } fprintf(stderr, "Generating 1-4 interactions: fudge = %g\n", fudge); - snew(pairs->param, pairs->nr); - for (int i = 0; (i < ntp); i++) + pairs->interactionTypes.clear(); + int i = 0; + std::array atomNumbers; + std::array forceParam = {NOTSET}; + for (const auto &type : nbs.interactionTypes) { - /* Copy param.a */ - pairs->param[i].a[0] = i / nnn; - pairs->param[i].a[1] = i % nnn; + /* Copy type.atoms */ + atomNumbers = {i / nnn, i % nnn }; /* Copy normal and FEP parameters and multiply by fudge factor */ - - - - for (int j = 0; (j < nrfp); j++) + gmx::ArrayRef existingParam = type.forceParam(); + GMX_RELEASE_ASSERT(2*nrfp <= MAXFORCEPARAM, "Can't have more parameters than half of maximum p arameter number"); + for (int j = 0; j < nrfp; j++) { /* If we are using sigma/epsilon values, only the epsilon values * should be scaled, but not sigma. @@ -127,13 +127,11 @@ static void gen_pairs(const InteractionTypeParameters &nbs, InteractionTypeParam scaling = fudge; } - 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]; + forceParam[j] = scaling*existingParam[j]; + forceParam[nrfp+j] = scaling*existingParam[j]; } + pairs->interactionTypes.emplace_back(InteractionType(atomNumbers, forceParam)); + i++; } } diff --git a/src/gromacs/gmxpreprocess/toppush.cpp b/src/gromacs/gmxpreprocess/toppush.cpp index 0d6dbcd7cc..c4d7b16420 100644 --- a/src/gromacs/gmxpreprocess/toppush.cpp +++ b/src/gromacs/gmxpreprocess/toppush.cpp @@ -72,16 +72,15 @@ void generate_nbparams(int comb, PreprocessingAtomTypes *atypes, warninp *wi) { - int k = -1; int nr, nrfp; real c, bi, bj, ci, cj, ci0, ci1, ci2, cj0, cj1, cj2; /* Lean mean shortcuts */ nr = atypes->size(); nrfp = NRFP(ftype); - snew(plist->param, nr*nr); - plist->nr = nr*nr; + plist->interactionTypes.clear(); + std::array forceParam = {NOTSET}; /* Fill the matrix with force parameters */ switch (ftype) { @@ -90,63 +89,66 @@ void generate_nbparams(int comb, { case eCOMB_GEOMETRIC: /* Gromos rules */ - for (int i = k = 0; (i < nr); i++) + for (int i = 0; (i < nr); i++) { - for (int j = 0; (j < nr); j++, k++) + for (int j = 0; (j < nr); j++) { for (int nf = 0; (nf < nrfp); nf++) { - ci = atypes->atomNonBondedParamFromAtomType(i, nf); - cj = atypes->atomNonBondedParamFromAtomType(j, nf); - c = std::sqrt(ci * cj); - plist->param[k].c[nf] = c; + ci = atypes->atomNonBondedParamFromAtomType(i, nf); + cj = atypes->atomNonBondedParamFromAtomType(j, nf); + c = std::sqrt(ci * cj); + forceParam[nf] = c; } + plist->interactionTypes.emplace_back(InteractionType({}, forceParam)); } } break; case eCOMB_ARITHMETIC: /* c0 and c1 are sigma and epsilon */ - for (int i = k = 0; (i < nr); i++) + for (int i = 0; (i < nr); i++) { - for (int j = 0; (j < nr); j++, k++) + for (int j = 0; (j < nr); j++) { 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; + forceParam[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. */ if (ci0 < 0 || cj0 < 0) { - plist->param[k].c[0] *= -1; + forceParam[0] *= -1; } - plist->param[k].c[1] = std::sqrt(ci1*cj1); + forceParam[1] = std::sqrt(ci1*cj1); + plist->interactionTypes.emplace_back(InteractionType({}, forceParam)); } } break; case eCOMB_GEOM_SIG_EPS: /* c0 and c1 are sigma and epsilon */ - for (int i = k = 0; (i < nr); i++) + for (int i = 0; (i < nr); i++) { - for (int j = 0; (j < nr); j++, k++) + for (int j = 0; (j < nr); j++) { 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)); + forceParam[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. */ if (ci0 < 0 || cj0 < 0) { - plist->param[k].c[0] *= -1; + forceParam[0] *= -1; } - plist->param[k].c[1] = std::sqrt(ci1*cj1); + forceParam[1] = std::sqrt(ci1*cj1); + plist->interactionTypes.emplace_back(InteractionType({}, forceParam)); } } @@ -155,17 +157,13 @@ void generate_nbparams(int comb, auto message = gmx::formatString("No such combination rule %d", comb); warning_error_and_exit(wi, message, FARGS); } - if (plist->nr != k) - { - gmx_incons("Topology processing, generate nb parameters"); - } break; case F_BHAM: /* Buckingham rules */ - for (int i = k = 0; (i < nr); i++) + for (int i = 0; (i < nr); i++) { - for (int j = 0; (j < nr); j++, k++) + for (int j = 0; (j < nr); j++) { ci0 = atypes->atomNonBondedParamFromAtomType(i, 0); cj0 = atypes->atomNonBondedParamFromAtomType(j, 0); @@ -173,16 +171,17 @@ void generate_nbparams(int comb, cj2 = atypes->atomNonBondedParamFromAtomType(j, 2); bi = atypes->atomNonBondedParamFromAtomType(i, 1); bj = atypes->atomNonBondedParamFromAtomType(j, 1); - plist->param[k].c[0] = std::sqrt(ci0 * cj0); + forceParam[0] = std::sqrt(ci0 * cj0); if ((bi == 0) || (bj == 0)) { - plist->param[k].c[1] = 0; + forceParam[1] = 0; } else { - plist->param[k].c[1] = 2.0/(1/bi+1/bj); + forceParam[1] = 2.0/(1/bi+1/bj); } - plist->param[k].c[2] = std::sqrt(ci2 * cj2); + forceParam[2] = std::sqrt(ci2 * cj2); + plist->interactionTypes.emplace_back(InteractionType({}, forceParam)); } } @@ -220,23 +219,22 @@ static void realloc_nb_params(PreprocessingAtomTypes *atypes, int copy_nbparams(t_nbparam **param, int ftype, InteractionTypeParameters *plist, int nr) { - int i, j, f; int nrfp, ncopy; nrfp = NRFP(ftype); ncopy = 0; - for (i = 0; i < nr; i++) + for (int i = 0; i < nr; i++) { - for (j = 0; j <= i; j++) + for (int j = 0; j <= i; j++) { GMX_RELEASE_ASSERT(param, "Must have valid parameters"); if (param[i][j].bSet) { - for (f = 0; f < nrfp; f++) + for (int f = 0; f < nrfp; f++) { - plist->param[nr*i+j].c[f] = param[i][j].c[f]; - plist->param[nr*j+i].c[f] = param[i][j].c[f]; + plist->interactionTypes[nr*i+j].setForceParameter(f, param[i][j].c[f]); + plist->interactionTypes[nr*j+i].setForceParameter(f, param[i][j].c[f]); } ncopy++; } @@ -282,7 +280,7 @@ void push_at (t_symtab *symtab, PreprocessingAtomTypes *at, gpp_bond_atomtype *b const char *entry; int ptype; } t_xlate; - t_xlate xl[eptNR] = { + t_xlate xl[eptNR] = { { "A", eptAtom }, { "N", eptNucleus }, { "S", eptShell }, @@ -290,20 +288,18 @@ void push_at (t_symtab *symtab, PreprocessingAtomTypes *at, gpp_bond_atomtype *b { "V", eptVSite }, }; - int nr, i, nfields, j, pt, nfp0 = -1; - int batype_nr, nread; - char type[STRLEN], btype[STRLEN], ptype[STRLEN]; - double m, q; - double c[MAXFORCEPARAM]; - char tmpfield[12][100]; /* Max 12 fields of width 100 */ - t_atom *atom; - t_param *param; - int atomnr; - bool have_atomic_number; - bool have_bonded_type; + int nr, nfields, j, pt, nfp0 = -1; + int batype_nr, nread; + char type[STRLEN], btype[STRLEN], ptype[STRLEN]; + double m, q; + double c[MAXFORCEPARAM]; + char tmpfield[12][100]; /* Max 12 fields of width 100 */ + t_atom *atom; + int atomnr; + bool have_atomic_number; + bool have_bonded_type; snew(atom, 1); - snew(param, 1); /* First assign input line to temporary array */ nfields = sscanf(line, "%s%s%s%s%s%s%s%s%s%s%s%s", @@ -506,6 +502,7 @@ void push_at (t_symtab *symtab, PreprocessingAtomTypes *at, gpp_bond_atomtype *b { c[j] = 0.0; } + std::array forceParam; if (strlen(type) == 1 && isdigit(type[0])) { @@ -539,11 +536,13 @@ void push_at (t_symtab *symtab, PreprocessingAtomTypes *at, gpp_bond_atomtype *b atom->q = q; atom->m = m; atom->ptype = pt; - for (i = 0; (i < MAXFORCEPARAM); i++) + for (int i = 0; i < MAXFORCEPARAM; i++) { - param->c[i] = c[i]; + forceParam[i] = c[i]; } + InteractionType interactionType({}, forceParam, ""); + if ((batype_nr = get_bond_atomtype_type(btype, bat)) == NOTSET) { add_bond_atomtype(bat, symtab, btype); @@ -554,14 +553,14 @@ void push_at (t_symtab *symtab, PreprocessingAtomTypes *at, gpp_bond_atomtype *b { auto message = gmx::formatString("Overriding atomtype %s", type); warning(wi, message); - if ((nr = at->setType(nr, symtab, *atom, type, param, batype_nr, + if ((nr = at->setType(nr, symtab, *atom, type, interactionType, batype_nr, atomnr)) == NOTSET) { auto message = gmx::formatString("Replacing atomtype %s failed", type); warning_error_and_exit(wi, message, FARGS); } } - else if ((at->addType(symtab, *atom, type, param, + else if ((at->addType(symtab, *atom, type, interactionType, batype_nr, atomnr)) == NOTSET) { auto message = gmx::formatString("Adding atomtype %s failed", type); @@ -573,7 +572,6 @@ void push_at (t_symtab *symtab, PreprocessingAtomTypes *at, gpp_bond_atomtype *b realloc_nb_params(at, nbparam, pair); } sfree(atom); - sfree(param); } //! Return whether the contents of \c a and \c b are the same, considering also reversed order. @@ -584,15 +582,15 @@ static bool equalEitherForwardOrBackward(gmx::ArrayRef a, gmx::ArrayRef std::equal(a.begin(), a.end(), b.rbegin())); } -static void push_bondtype(InteractionTypeParameters * bt, - const t_param * b, - int nral, - int ftype, - bool bAllowRepeat, - const char * line, - warninp *wi) +static void push_bondtype(InteractionTypeParameters *bt, + const InteractionType &b, + int nral, + int ftype, + bool bAllowRepeat, + const char *line, + warninp *wi) { - int nr = bt->nr; + int nr = bt->size(); int nrfp = NRFP(ftype); /* If bAllowRepeat is TRUE, we allow multiple entries as long as they @@ -609,9 +607,11 @@ static void push_bondtype(InteractionTypeParameters * bt, if (bAllowRepeat && nr > 1) { isContinuationOfBlock = true; + gmx::ArrayRef newParAtom = b.atoms(); + gmx::ArrayRef sysParAtom = bt->interactionTypes[nr - 2].atoms(); for (int j = 0; j < nral; j++) { - if (b->a[j] != bt->param[nr - 2].a[j]) + if (newParAtom[j] != sysParAtom[j]) { isContinuationOfBlock = false; } @@ -626,12 +626,14 @@ static void push_bondtype(InteractionTypeParameters * bt, bool haveErrored = false; for (int i = 0; (i < nr); i++) { - gmx::ArrayRef bParams(b->a, b->a + nral); - gmx::ArrayRef testParams(bt->param[i].a, bt->param[i].a + nral); + gmx::ArrayRef bParams = b.atoms(); + gmx::ArrayRef testParams = bt->interactionTypes[i].atoms(); + GMX_RELEASE_ASSERT(bParams.size() == testParams.size(), "Number of atoms needs to be the same between parameters"); if (equalEitherForwardOrBackward(bParams, testParams)) { GMX_ASSERT(nrfp <= MAXFORCEPARAM, "This is ensured in other places, but we need this assert to keep the clang analyzer happy"); - const bool identicalParameters = std::equal(bt->param[i].c, bt->param[i].c + nrfp, b->c); + const bool identicalParameters = std::equal(bt->interactionTypes[i].forceParam().begin(), + bt->interactionTypes[i].forceParam().begin() + nrfp, b.forceParam().begin()); if (!bAllowRepeat || identicalParameters) { @@ -666,9 +668,10 @@ static void push_bondtype(InteractionTypeParameters * bt, warning(wi, message); fprintf(stderr, " old: "); - for (int j = 0; (j < nrfp); j++) + gmx::ArrayRef forceParam = bt->interactionTypes[i].forceParam(); + for (int j = 0; j < nrfp; j++) { - fprintf(stderr, " %g", bt->param[i].c[j]); + fprintf(stderr, " %g", forceParam[j]); } fprintf(stderr, " \n new: %s\n\n", line); @@ -681,9 +684,10 @@ static void push_bondtype(InteractionTypeParameters * bt, /* Overwrite the parameters with the latest ones */ // TODO considering improving the following code by replacing with: // std::copy(b->c, b->c + nrfp, bt->param[i].c); - for (int j = 0; (j < nrfp); j++) + gmx::ArrayRef forceParam = b.forceParam(); + for (int j = 0; j < nrfp; j++) { - bt->param[i].c[j] = b->c[j]; + bt->interactionTypes[i].setForceParameter(j, forceParam[j]); } } } @@ -691,13 +695,11 @@ static void push_bondtype(InteractionTypeParameters * bt, if (addBondType) { - /* alloc */ - pr_alloc (2, bt); - /* fill the arrays up and down */ - memcpy(bt->param[bt->nr].c, b->c, sizeof(b->c)); - memcpy(bt->param[bt->nr].a, b->a, sizeof(b->a)); - memcpy(bt->param[bt->nr+1].c, b->c, sizeof(b->c)); + bt->interactionTypes.emplace_back( + InteractionType(b.atoms(), b.forceParam(), b.interactionTypeName())); + /* need to store force values because they might change below */ + std::vector forceParam(b.forceParam().begin(), b.forceParam().end()); /* The definitions of linear angles depend on the order of atoms, * that means that for atoms i-j-k, with certain parameter a, the @@ -705,16 +707,16 @@ static void push_bondtype(InteractionTypeParameters * bt, */ if (ftype == F_LINEAR_ANGLES) { - bt->param[bt->nr+1].c[0] = 1-bt->param[bt->nr+1].c[0]; - bt->param[bt->nr+1].c[2] = 1-bt->param[bt->nr+1].c[2]; + forceParam[0] = 1- forceParam[0]; + forceParam[2] = 1- forceParam[2]; } - - for (int j = 0; (j < nral); j++) + std::vector atoms; + gmx::ArrayRef oldAtoms = b.atoms(); + for (auto oldAtom = oldAtoms.rbegin(); oldAtom != oldAtoms.rend(); oldAtom++) { - bt->param[bt->nr+1].a[j] = b->a[nral-1-j]; + atoms.emplace_back(*oldAtom); } - - bt->nr += 2; + bt->interactionTypes.emplace_back(InteractionType(atoms, forceParam, b.interactionTypeName())); } } @@ -726,7 +728,7 @@ void push_bt(Directive d, char *line, warninp *wi) { - const char *formal[MAXATOMLIST+1] = { + const char *formal[MAXATOMLIST+1] = { "%s", "%s%s", "%s%s%s", @@ -735,7 +737,7 @@ void push_bt(Directive d, "%s%s%s%s%s%s", "%s%s%s%s%s%s%s" }; - const char *formnl[MAXATOMLIST+1] = { + const char *formnl[MAXATOMLIST+1] = { "%*s", "%*s%*s", "%*s%*s%*s", @@ -744,13 +746,12 @@ void push_bt(Directive d, "%*s%*s%*s%*s%*s%*s", "%*s%*s%*s%*s%*s%*s%*s" }; - const char *formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf"; - int i, ft, ftype, nn, nrfp, nrfpA; - char f1[STRLEN]; - char alc[MAXATOMLIST+1][20]; + const char *formlf = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf"; + int i, ft, ftype, nn, nrfp, nrfpA; + char f1[STRLEN]; + char alc[MAXATOMLIST+1][20]; /* One force parameter more, so we can check if we read too many */ - double c[MAXFORCEPARAM+1]; - t_param p; + double c[MAXFORCEPARAM+1]; if ((bat && at) || (!bat && !at)) { @@ -800,24 +801,28 @@ void push_bt(Directive d, } } } - for (i = 0; (i < nral); i++) + std::vector atoms; + std::array forceParam; + for (int i = 0; (i < nral); i++) { - if ((at != nullptr) && ((p.a[i] = at->atomTypeFromName(alc[i])) == NOTSET)) + int atomNumber; + if ((at != nullptr) && ((atomNumber = at->atomTypeFromName(alc[i])) == NOTSET)) { auto message = gmx::formatString("Unknown atomtype %s\n", alc[i]); warning_error_and_exit(wi, message, FARGS); } - else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET)) + else if (bat && ((atomNumber = get_bond_atomtype_type(alc[i], bat)) == NOTSET)) { auto message = gmx::formatString("Unknown bond_atomtype %s\n", alc[i]); warning_error_and_exit(wi, message, FARGS); } + atoms.emplace_back(atomNumber); } - for (i = 0; (i < nrfp); i++) + for (int i = 0; (i < nrfp); i++) { - p.c[i] = c[i]; + forceParam[i] = c[i]; } - push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE, line, wi); + push_bondtype (&(bt[ftype]), InteractionType(atoms, forceParam), nral, ftype, FALSE, line, wi); } @@ -825,7 +830,7 @@ void push_dihedraltype(Directive d, gmx::ArrayRef bt, gpp_bond_atomtype *bat, char *line, warninp *wi) { - const char *formal[MAXATOMLIST+1] = { + const char *formal[MAXATOMLIST+1] = { "%s", "%s%s", "%s%s%s", @@ -834,7 +839,7 @@ void push_dihedraltype(Directive d, gmx::ArrayRef bt, "%s%s%s%s%s%s", "%s%s%s%s%s%s%s" }; - const char *formnl[MAXATOMLIST+1] = { + const char *formnl[MAXATOMLIST+1] = { "%*s", "%*s%*s", "%*s%*s%*s", @@ -843,7 +848,7 @@ void push_dihedraltype(Directive d, gmx::ArrayRef bt, "%*s%*s%*s%*s%*s%*s", "%*s%*s%*s%*s%*s%*s%*s" }; - const char *formlf[MAXFORCEPARAM] = { + const char *formlf[MAXFORCEPARAM] = { "%lf", "%lf%lf", "%lf%lf%lf", @@ -857,12 +862,11 @@ void push_dihedraltype(Directive d, gmx::ArrayRef bt, "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf", "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf", }; - int i, ft, ftype, nn, nrfp, nrfpA, nral; - char f1[STRLEN]; - char alc[MAXATOMLIST+1][20]; - double c[MAXFORCEPARAM]; - t_param p; - bool bAllowRepeat; + int i, ft, ftype, nn, nrfp, nrfpA, nral; + char f1[STRLEN]; + char alc[MAXATOMLIST+1][20]; + double c[MAXFORCEPARAM]; + bool bAllowRepeat; /* This routine accepts dihedraltypes defined from either 2 or 4 atoms. * @@ -964,29 +968,33 @@ void push_dihedraltype(Directive d, gmx::ArrayRef bt, } } - for (i = 0; (i < 4); i++) + std::vector atoms; + std::array forceParam; + for (int i = 0; (i < 4); i++) { if (!strcmp(alc[i], "X")) { - p.a[i] = -1; + atoms.emplace_back(-1); } else { - if ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET) + int atomNumber; + if ((atomNumber = get_bond_atomtype_type(alc[i], bat)) == NOTSET) { auto message = gmx::formatString("Unknown bond_atomtype %s", alc[i]); warning_error_and_exit(wi, message, FARGS); } + atoms.emplace_back(atomNumber); } } - for (i = 0; (i < nrfp); i++) + for (int i = 0; (i < nrfp); i++) { - p.c[i] = c[i]; + forceParam[i] = c[i]; } /* Always use 4 atoms here, since we created two wildcard atoms * if there wasn't of them 4 already. */ - push_bondtype (&(bt[ftype]), &p, 4, ftype, bAllowRepeat, line, wi); + push_bondtype (&(bt[ftype]), InteractionType(atoms, forceParam), 4, ftype, bAllowRepeat, line, wi); } @@ -1125,13 +1133,12 @@ push_cmaptype(Directive d, char *line, warninp *wi) { - const char *formal = "%s%s%s%s%s%s%s%s%n"; + const char *formal = "%s%s%s%s%s%s%s%s%n"; - 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]; - t_param p; + 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]; /* Keep the compiler happy */ read_cmap = 0; @@ -1221,18 +1228,21 @@ push_cmaptype(Directive d, bt[F_CMAP].cmapAngles++; /* Since we are incrementing here, we need to subtract later, see (*****) */ nct = (nral+1) * bt[F_CMAP].cmapAngles; + std::vector atoms; for (int i = 0; (i < nral); i++) { - if (at && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET)) + int atomNumber; + if (at && ((atomNumber = get_bond_atomtype_type(alc[i], bat)) == NOTSET)) { auto message = gmx::formatString("Unknown atomtype %s\n", alc[i]); warning_error(wi, message); } - else if (bat && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET)) + else if (bat && ((atomNumber = get_bond_atomtype_type(alc[i], bat)) == NOTSET)) { auto message = gmx::formatString("Unknown bond_atomtype %s\n", alc[i]); warning_error(wi, message); } + atoms.emplace_back(atomNumber); /* Assign a grid number to each cmap_type */ bt[F_CMAP].cmapAtomTypes.emplace_back(get_bond_atomtype_type(alc[i], bat)); @@ -1248,14 +1258,10 @@ push_cmaptype(Directive d, warning_error(wi, message); } - /* Is this correct?? */ - for (int i = 0; i < MAXFORCEPARAM; i++) - { - p.c[i] = NOTSET; - } + std::array forceParam = {NOTSET}; /* Push the bond to the bondlist */ - push_bondtype (&(bt[ftype]), &p, nral, ftype, FALSE, line, wi); + push_bondtype (&(bt[ftype]), InteractionType(atoms, forceParam), nral, ftype, FALSE, line, wi); } @@ -1472,16 +1478,49 @@ void push_molt(t_symtab *symtab, mol->back().excl_set = false; } +static bool findIfAllNBAtomsMatch(gmx::ArrayRef atomsFromParameterArray, + gmx::ArrayRef atomsFromCurrentParameter, + const t_atoms *at, + bool bB) +{ + if (atomsFromParameterArray.size() != atomsFromCurrentParameter.size()) + { + return false; + } + else if (bB) + { + for (int i = 0; i < atomsFromCurrentParameter.ssize(); i++) + { + if (at->atom[atomsFromCurrentParameter[i]].typeB != atomsFromParameterArray[i]) + { + return false; + } + } + return true; + } + else + { + for (int i = 0; i < atomsFromCurrentParameter.ssize(); i++) + { + if (at->atom[atomsFromCurrentParameter[i]].type != atomsFromParameterArray[i]) + { + return false; + } + } + return true; + } +} + static bool default_nb_params(int ftype, gmx::ArrayRef bt, t_atoms *at, - t_param *p, int c_start, bool bB, bool bGenPairs) + InteractionType *p, int c_start, bool bB, bool bGenPairs) { - int i, j, ti, tj, ntype; - bool bFound; - t_param *pi = nullptr; - int nr = bt[ftype].nr; - int nral = NRAL(ftype); - int nrfp = interaction_function[ftype].nrfpA; - int nrfpB = interaction_function[ftype].nrfpB; + int ti, tj, ntype; + bool bFound; + InteractionType *pi = nullptr; + int nr = bt[ftype].size(); + int nral = NRAL(ftype); + int nrfp = interaction_function[ftype].nrfpA; + int nrfpB = interaction_function[ftype].nrfpB; if ((!bB && nrfp == 0) || (bB && nrfpB == 0)) { @@ -1498,70 +1537,68 @@ static bool default_nb_params(int ftype, gmx::ArrayRefatom[p->a[0]].typeB; - tj = at->atom[p->a[1]].typeB; + ti = at->atom[p->ai()].typeB; + tj = at->atom[p->aj()].typeB; + } + else + { + ti = at->atom[p->ai()].type; + tj = at->atom[p->aj()].type; + } + pi = &(bt[ftype].interactionTypes[ntype*ti+tj]); + if (pi->atoms().ssize() < nral) + { + /* not initialized yet with atom names */ + bFound = false; } else { - ti = at->atom[p->a[0]].type; - tj = at->atom[p->a[1]].type; + bFound = ((ti == pi->ai()) && (tj == pi->aj())); } - pi = &(bt[ftype].param[ntype*ti+tj]); - bFound = ((ti == pi->a[0]) && (tj == pi->a[1])); } + gmx::ArrayRef paramAtoms = p->atoms(); /* Search explicitly if we didnt find it */ if (!bFound) { - for (i = 0; ((i < nr) && !bFound); i++) + auto foundParameter = std::find_if(bt[ftype].interactionTypes.begin(), + bt[ftype].interactionTypes.end(), + [¶mAtoms, &at, &bB](const auto ¶m) + { return findIfAllNBAtomsMatch(param.atoms(), paramAtoms, at, bB); }); + if (foundParameter != bt[ftype].interactionTypes.end()) { - pi = &(bt[ftype].param[i]); - if (bB) - { - for (j = 0; ((j < nral) && - (at->atom[p->a[j]].typeB == pi->a[j])); j++) - { - ; - } - } - else - { - for (j = 0; ((j < nral) && - (at->atom[p->a[j]].type == pi->a[j])); j++) - { - ; - } - } - bFound = (j == nral); + bFound = true; + pi = &(*foundParameter); } } if (bFound) { + gmx::ArrayRef forceParam = pi->forceParam(); if (bB) { if (nrfp+nrfpB > MAXFORCEPARAM) { gmx_incons("Too many force parameters"); } - for (j = c_start; (j < nrfpB); j++) + for (int j = c_start; j < nrfpB; j++) { - p->c[nrfp+j] = pi->c[j]; + p->setForceParameter(nrfp+j, forceParam[j]); } } else { - for (j = c_start; (j < nrfp); j++) + for (int j = c_start; j < nrfp; j++) { - p->c[j] = pi->c[j]; + p->setForceParameter(j, forceParam[j]); } } } else { - for (j = c_start; (j < nrfp); j++) + for (int j = c_start; j < nrfp; j++) { - p->c[j] = 0.0; + p->setForceParameter(j, 0.0); } } return bFound; @@ -1569,7 +1606,7 @@ static bool default_nb_params(int ftype, gmx::ArrayRef bondtype, t_atoms *at, PreprocessingAtomTypes *atypes, - t_param *p, bool bB, + InteractionType *p, bool bB, int *cmap_type, int *nparam_def, warninp *wi) { @@ -1590,11 +1627,11 @@ static bool default_cmap_params(gmx::ArrayRef bondtyp else { if ( - (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])) + (atypes->bondAtomTypeFromAtomType(at->atom[p->ai()].type) == bondtype[F_CMAP].cmapAtomTypes[i]) && + (atypes->bondAtomTypeFromAtomType(at->atom[p->aj()].type) == bondtype[F_CMAP].cmapAtomTypes[i+1]) && + (atypes->bondAtomTypeFromAtomType(at->atom[p->ak()].type) == bondtype[F_CMAP].cmapAtomTypes[i+2]) && + (atypes->bondAtomTypeFromAtomType(at->atom[p->al()].type) == bondtype[F_CMAP].cmapAtomTypes[i+3]) && + (atypes->bondAtomTypeFromAtomType(at->atom[p->am()].type) == bondtype[F_CMAP].cmapAtomTypes[i+4])) { /* Found cmap torsion */ bFound = true; @@ -1608,7 +1645,7 @@ static bool default_cmap_params(gmx::ArrayRef bondtyp if (!bFound) { auto message = gmx::formatString("Unknown cmap torsion between atoms %d %d %d %d %d", - p->a[0]+1, p->a[1]+1, p->a[2]+1, p->a[3]+1, p->a[4]+1); + p->ai()+1, p->aj()+1, p->ak()+1, p->al()+1, p->am()+1); warning_error_and_exit(wi, message, FARGS); } @@ -1621,20 +1658,20 @@ static bool default_cmap_params(gmx::ArrayRef bondtyp /* Returns the number of exact atom type matches, i.e. non wild-card matches, * returns -1 when there are no matches at all. */ -static int natom_match(t_param *pi, +static int natom_match(const InteractionType &pi, int type_i, int type_j, int type_k, int type_l, const PreprocessingAtomTypes* atypes) { - 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())) + 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) + - (pi->aj() == -1 ? 0 : 1) + - (pi->ak() == -1 ? 0 : 1) + - (pi->al() == -1 ? 0 : 1); + (pi.ai() == -1 ? 0 : 1) + + (pi.aj() == -1 ? 0 : 1) + + (pi.ak() == -1 ? 0 : 1) + + (pi.al() == -1 ? 0 : 1); } else { @@ -1642,65 +1679,106 @@ static int natom_match(t_param *pi, } } -static bool defaultInteractionTypeParameters(int ftype, gmx::ArrayRef bt, - t_atoms *at, PreprocessingAtomTypes *atypes, - t_param *p, bool bB, - t_param **param_def, - int *nparam_def) +static int findNumberOfDihedralAtomMatches(const InteractionType ¤tParamFromParameterArray, + const InteractionType ¶meterToAdd, + const t_atoms *at, + const PreprocessingAtomTypes *atypes, + bool bB) { - int nparam_found; - bool bFound, bSame; - t_param *pi = nullptr; - t_param *pj = nullptr; - int nr = bt[ftype].nr; - int nral = NRAL(ftype); - int nrfpA = interaction_function[ftype].nrfpA; - int nrfpB = interaction_function[ftype].nrfpB; + if (bB) + { + return natom_match(currentParamFromParameterArray, + at->atom[parameterToAdd.ai()].typeB, + at->atom[parameterToAdd.aj()].typeB, + at->atom[parameterToAdd.ak()].typeB, + at->atom[parameterToAdd.al()].typeB, atypes); + } + else + { + return natom_match(currentParamFromParameterArray, + at->atom[parameterToAdd.ai()].type, + at->atom[parameterToAdd.aj()].type, + at->atom[parameterToAdd.ak()].type, + at->atom[parameterToAdd.al()].type, atypes); + } +} + +static bool findIfAllParameterAtomsMatch(gmx::ArrayRef atomsFromParameterArray, + gmx::ArrayRef atomsFromCurrentParameter, + const t_atoms *at, + const PreprocessingAtomTypes *atypes, + bool bB) +{ + if (atomsFromParameterArray.size() != atomsFromCurrentParameter.size()) + { + return false; + } + else if (bB) + { + for (int i = 0; i < atomsFromCurrentParameter.ssize(); i++) + { + if (atypes->bondAtomTypeFromAtomType( + at->atom[atomsFromCurrentParameter[i]].typeB) != atomsFromParameterArray[i]) + { + return false; + } + } + return true; + } + else + { + for (int i = 0; i < atomsFromCurrentParameter.ssize(); i++) + { + if (atypes->bondAtomTypeFromAtomType( + at->atom[atomsFromCurrentParameter[i]].type) != atomsFromParameterArray[i]) + { + return false; + } + } + return true; + } +} + +static std::vector::iterator +defaultInteractionTypeParameters(int ftype, gmx::ArrayRef bt, + t_atoms *at, PreprocessingAtomTypes *atypes, + const InteractionType &p, bool bB, + int *nparam_def) +{ + int nparam_found; + int nrfpA = interaction_function[ftype].nrfpA; + int nrfpB = interaction_function[ftype].nrfpB; if ((!bB && nrfpA == 0) || (bB && nrfpB == 0)) { - return TRUE; + return bt[ftype].interactionTypes.end(); } - bFound = FALSE; nparam_found = 0; if (ftype == F_PDIHS || ftype == F_RBDIHS || ftype == F_IDIHS || ftype == F_PIDIHS) { int nmatch_max = -1; - int i = -1; - int t; /* For dihedrals we allow wildcards. We choose the first type * that has the most real matches, i.e. non-wildcard matches. */ - for (t = 0; ((t < nr) && nmatch_max < 4); t++) - { - int nmatch; - t_param *pt; - - pt = &(bt[ftype].param[t]); - if (bB) - { - 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, atypes); - } - if (nmatch > nmatch_max) + auto prevPos = bt[ftype].interactionTypes.end(); + auto pos = bt[ftype].interactionTypes.begin(); + while (pos != bt[ftype].interactionTypes.end() && nmatch_max < 4) + { + pos = std::find_if(bt[ftype].interactionTypes.begin(), bt[ftype].interactionTypes.end(), + [&p, &at, &atypes, &bB, &nmatch_max](const auto ¶m) + { return (findNumberOfDihedralAtomMatches(param, p, at, atypes, bB) > nmatch_max); }); + if (pos != bt[ftype].interactionTypes.end()) { - nmatch_max = nmatch; - i = t; - bFound = TRUE; + prevPos = pos; + nmatch_max = findNumberOfDihedralAtomMatches(*pos, p, at, atypes, bB); } } - if (bFound) + if (prevPos != bt[ftype].interactionTypes.end()) { - int j; - - pi = &(bt[ftype].param[i]); nparam_found++; /* Find additional matches for this dihedral - necessary @@ -1708,12 +1786,11 @@ static bool defaultInteractionTypeParameters(int ftype, gmx::ArrayRefai() == pj->ai() && pi->aj() == pj->aj() && pi->ak() == pj->ak() && pi->al() == pj->al()); + bSame = (prevPos->ai() == nextPos->ai() && prevPos->aj() == nextPos->aj() && prevPos->ak() == nextPos->ak() && prevPos->al() == nextPos->al()); if (bSame) { nparam_found++; @@ -1721,42 +1798,23 @@ static bool defaultInteractionTypeParameters(int ftype, gmx::ArrayRefbondAtomTypeFromAtomType(at->atom[p->a[j]].typeB) == pi->a[j])); j++) - { - ; - } - } - else - { - for (j = 0; ((j < nral) && - (atypes->bondAtomTypeFromAtomType(at->atom[p->a[j]].type) == pi->a[j])); j++) - { - ; - } - } - bFound = (j == nral); - } - if (bFound) + gmx::ArrayRef atomParam = p.atoms(); + auto found = std::find_if(bt[ftype].interactionTypes.begin(), + bt[ftype].interactionTypes.end(), + [&atomParam, &at, &atypes, &bB](const auto ¶m) + { return findIfAllParameterAtomsMatch(param.atoms(), atomParam, at, atypes, bB); }); + if (found != bt[ftype].interactionTypes.end()) { nparam_found = 1; } + *nparam_def = nparam_found; + return found; } - - *param_def = pi; - *nparam_def = nparam_found; - - return bFound; } @@ -1768,7 +1826,7 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, bool bZero, bool *bWarn_copy_A_B, warninp *wi) { - const char *aaformat[MAXATOMLIST] = { + const char *aaformat[MAXATOMLIST] = { "%d%d", "%d%d%d", "%d%d%d%d", @@ -1776,7 +1834,7 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, "%d%d%d%d%d%d", "%d%d%d%d%d%d%d" }; - const char *asformat[MAXATOMLIST] = { + const char *asformat[MAXATOMLIST] = { "%*s%*s", "%*s%*s%*s", "%*s%*s%*s%*s", @@ -1784,21 +1842,20 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, "%*s%*s%*s%*s%*s%*s", "%*s%*s%*s%*s%*s%*s%*s" }; - const char *ccformat = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf"; - int nr, i, j, nral, nral_fmt, nread, ftype; - char format[STRLEN]; + const char *ccformat = "%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf%lf"; + int nral, nral_fmt, nread, ftype; + char format[STRLEN]; /* One force parameter more, so we can check if we read too many */ - double cc[MAXFORCEPARAM+1]; - int aa[MAXATOMLIST+1]; - t_param param, *param_defA, *param_defB; - bool bFoundA = FALSE, bFoundB = FALSE, bDef, bPert, bSwapParity = FALSE; - int nparam_defA, nparam_defB; + double cc[MAXFORCEPARAM+1]; + int aa[MAXATOMLIST+1]; + bool bFoundA = FALSE, bFoundB = FALSE, bDef, bSwapParity = FALSE; + int nparam_defA, nparam_defB; nparam_defA = nparam_defB = 0; ftype = ifunc_index(d, 1); nral = NRAL(ftype); - for (j = 0; j < MAXATOMLIST; j++) + for (int j = 0; j < nral; j++) { aa[j] = NOTSET; } @@ -1858,7 +1915,7 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, /* Check for double atoms and atoms out of bounds */ - for (i = 0; (i < nral); i++) + for (int i = 0; (i < nral); i++) { if (aa[i] < 1 || aa[i] > at->nr) { @@ -1870,7 +1927,7 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, aa[i], dir2str(d), at->nr, dir2str(d), dir2str(d)); warning_error_and_exit(wi, message, FARGS); } - for (j = i+1; (j < nral); j++) + for (int j = i+1; (j < nral); j++) { GMX_ASSERT(j < MAXATOMLIST + 1, "Values from nral=NRAL() will satisfy this, we assert to keep gcc 4 happy"); if (aa[i] == aa[j]) @@ -1893,39 +1950,65 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, } /* default force parameters */ - for (j = 0; (j < MAXATOMLIST); j++) + std::vector atoms; + for (int j = 0; (j < nral); j++) { - param.a[j] = aa[j]-1; - } - for (j = 0; (j < MAXFORCEPARAM); j++) - { - param.c[j] = 0.0; + atoms.emplace_back(aa[j]-1); } + /* need to have an empty but initialized param array for some reason */ + std::array forceParam = {0.0}; /* Get force params for normal and free energy perturbation * studies, as determined by types! */ + InteractionType param(atoms, forceParam, ""); + std::vector::iterator foundAParameter = bondtype[ftype].interactionTypes.end(); + std::vector::iterator foundBParameter = bondtype[ftype].interactionTypes.end(); if (bBonded) { - bFoundA = defaultInteractionTypeParameters(ftype, bondtype, at, atypes, ¶m, FALSE, ¶m_defA, &nparam_defA); - if (bFoundA) + foundAParameter = defaultInteractionTypeParameters(ftype, + bondtype, + at, + atypes, + param, + FALSE, + &nparam_defA); + if (foundAParameter != bondtype[ftype].interactionTypes.end()) { /* Copy the A-state and B-state default parameters. */ GMX_ASSERT(NRFPA(ftype)+NRFPB(ftype) <= MAXFORCEPARAM, "Bonded interactions may have at most 12 parameters"); - for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++) + gmx::ArrayRef defaultParam = foundAParameter->forceParam(); + for (int j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++) { - param.c[j] = param_defA->c[j]; + param.setForceParameter(j, defaultParam[j]); } + bFoundA = true; + } + else if (NRFPA(ftype) == 0) + { + bFoundA = true; } - bFoundB = defaultInteractionTypeParameters(ftype, bondtype, at, atypes, ¶m, TRUE, ¶m_defB, &nparam_defB); - if (bFoundB) + foundBParameter = defaultInteractionTypeParameters(ftype, + bondtype, + at, + atypes, + param, + TRUE, + &nparam_defB); + if (foundBParameter != bondtype[ftype].interactionTypes.end()) { /* Copy only the B-state default parameters */ - for (j = NRFPA(ftype); (j < NRFP(ftype)); j++) + gmx::ArrayRef defaultParam = foundBParameter->forceParam(); + for (int j = NRFPA(ftype); (j < NRFP(ftype)); j++) { - param.c[j] = param_defB->c[j]; + param.setForceParameter(j, defaultParam[j]); } + bFoundB = true; + } + else if (NRFPB(ftype) == 0) + { + bFoundB = true; } } else if (ftype == F_LJ14) @@ -1935,10 +2018,10 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, } else if (ftype == F_LJC14_Q) { - param.c[0] = fudgeQQ; /* Fill in the A-state charges as default parameters */ - param.c[1] = at->atom[param.a[0]].q; - param.c[2] = at->atom[param.a[1]].q; + param.setForceParameter(0, fudgeQQ); + param.setForceParameter(1, at->atom[param.ai()].q); + param.setForceParameter(2, at->atom[param.aj()].q); /* The default LJ parameters are the standard 1-4 parameters */ bFoundA = default_nb_params(F_LJ14, bondtype, at, ¶m, 3, FALSE, bGenPairs); bFoundB = TRUE; @@ -1967,10 +2050,11 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, if ((nread == NRFPA(ftype)) && (NRFPB(ftype) != 0)) { /* We only have to issue a warning if these atoms are perturbed! */ - bPert = FALSE; - for (j = 0; (j < nral); j++) + bool bPert = false; + gmx::ArrayRef paramAtoms = param.atoms(); + for (int j = 0; (j < nral); j++) { - bPert = bPert || PERTURBED(at->atom[param.a[j]]); + bPert = bPert || PERTURBED(at->atom[paramAtoms[j]]); } if (bPert && *bWarn_copy_A_B) @@ -1986,7 +2070,7 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, /* The B-state parameters correspond to the first nrfpB * A-state parameters. */ - for (j = 0; (j < NRFPB(ftype)); j++) + for (int j = 0; (j < NRFPB(ftype)); j++) { cc[nread++] = cc[j]; } @@ -2008,11 +2092,10 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, warning_error_and_exit(wi, message, FARGS); } - for (j = 0; (j < nread); j++) + for (int j = 0; (j < nread); j++) { - param.c[j] = cc[j]; + param.setForceParameter(j, cc[j]); } - /* Check whether we have to use the defaults */ if (nread == NRFP(ftype)) { @@ -2031,7 +2114,7 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, /* When we have multiple terms it would be very dangerous to allow perturbations to a different atom type! */ if (ftype == F_PDIHS) { - if ((nparam_defA != nparam_defB) || ((nparam_defA > 1 || nparam_defB > 1) && (param_defA != param_defB))) + if ((nparam_defA != nparam_defB) || ((nparam_defA > 1 || nparam_defB > 1) && (foundAParameter != foundBParameter))) { auto message = gmx::formatString("Cannot automatically perturb a torsion with multiple terms to different form.\n" @@ -2053,15 +2136,14 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, { if (interaction_function[ftype].flags & IF_VSITE) { - /* set them to NOTSET, will be calculated later */ - for (j = 0; (j < MAXFORCEPARAM); j++) + for (int j = 0; j < MAXFORCEPARAM; j++) { - param.c[j] = NOTSET; + param.setForceParameter(j, NOTSET); } - if (bSwapParity) { - param.c1() = -1; /* flag to swap parity of vsite construction */ + /* flag to swap parity of vsi te construction */ + param.setForceParameter(1, -1); } } else @@ -2085,10 +2167,10 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, switch (ftype) { case F_VSITE3FAD: - param.c0() = 360 - param.c0(); + param.setForceParameter(0, 360 - param.c0()); break; case F_VSITE3OUT: - param.c2() = -param.c2(); + param.setForceParameter(2, -param.c2()); break; } } @@ -2096,10 +2178,11 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, if (!bFoundB) { /* We only have to issue a warning if these atoms are perturbed! */ - bPert = FALSE; - for (j = 0; (j < nral); j++) + bool bPert = false; + gmx::ArrayRef paramAtoms = param.atoms(); + for (int j = 0; (j < nral); j++) { - bPert = bPert || PERTURBED(at->atom[param.a[j]]); + bPert = bPert || PERTURBED(at->atom[paramAtoms[j]]); } if (bPert) @@ -2112,30 +2195,32 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, } } + gmx::ArrayRef paramValue = param.forceParam(); if ((ftype == F_PDIHS || ftype == F_ANGRES || ftype == F_ANGRESZ) - && param.c[5] != param.c[2]) + && paramValue[5] != paramValue[2]) { auto message = gmx::formatString("%s multiplicity can not be perturbed %f!=%f", interaction_function[ftype].longname, - param.c[2], param.c[5]); + paramValue[2], paramValue[5]); warning_error_and_exit(wi, message, FARGS); } - if (IS_TABULATED(ftype) && param.c[2] != param.c[0]) + if (IS_TABULATED(ftype) && param.c0() != param.c2()) { auto message = gmx::formatString("%s table number can not be perturbed %d!=%d", interaction_function[ftype].longname, - gmx::roundToInt(param.c[0]), gmx::roundToInt(param.c[2])); + gmx::roundToInt(param.c0()), gmx::roundToInt(param.c0())); warning_error_and_exit(wi, message, FARGS); } /* Dont add R-B dihedrals where all parameters are zero (no interaction) */ if (ftype == F_RBDIHS) { - nr = 0; - for (i = 0; i < NRFP(ftype); i++) + + int nr = 0; + for (int i = 0; i < NRFP(ftype); i++) { - if (param.c[i] != 0) + if (paramValue[i] != 0.0) { nr++; } @@ -2147,7 +2232,7 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, } /* Put the values in the appropriate arrays */ - add_param_to_list (&bond[ftype], ¶m); + add_param_to_list (&bond[ftype], param); /* Push additional torsions from FF for ftype==9 if we have them. * We have already checked that the A/B states do not differ in this case, @@ -2156,16 +2241,17 @@ void push_bond(Directive d, gmx::ArrayRef bondtype, */ if (bDef && ftype == F_PDIHS) { - for (i = 1; i < nparam_defA; i++) + for (int i = 1; i < nparam_defA; i++) { /* Advance pointer! */ - param_defA += 2; - for (j = 0; (j < NRFPA(ftype)+NRFPB(ftype)); j++) + foundAParameter += 2; + gmx::ArrayRef forceParam = foundAParameter->forceParam(); + for (int j = 0; j < (NRFPA(ftype)+NRFPB(ftype)); j++) { - param.c[j] = param_defA->c[j]; + param.setForceParameter(j, forceParam[j]); } /* And push the next term for this torsion */ - add_param_to_list (&bond[ftype], ¶m); + add_param_to_list (&bond[ftype], param); } } } @@ -2186,11 +2272,10 @@ void push_cmap(Directive d, gmx::ArrayRef bondtype, "%d%d%d%d%d%d%d" }; - int i, j, ftype, nral, nread, ncmap_params; - int cmap_type; - int aa[MAXATOMLIST+1]; - bool bFound; - t_param param; + int ftype, nral, nread, ncmap_params; + int cmap_type; + int aa[MAXATOMLIST+1]; + bool bFound; ftype = ifunc_index(d, 1); nral = NRAL(ftype); @@ -2210,7 +2295,7 @@ void push_cmap(Directive d, gmx::ArrayRef bondtype, } /* Check for double atoms and atoms out of bounds */ - for (i = 0; i < nral; i++) + for (int i = 0; i < nral; i++) { if (aa[i] < 1 || aa[i] > at->nr) { @@ -2223,7 +2308,7 @@ void push_cmap(Directive d, gmx::ArrayRef bondtype, warning_error_and_exit(wi, message, FARGS); } - for (j = i+1; (j < nral); j++) + for (int j = i+1; (j < nral); j++) { if (aa[i] == aa[j]) { @@ -2234,15 +2319,13 @@ void push_cmap(Directive d, gmx::ArrayRef bondtype, } /* default force parameters */ - for (j = 0; (j < MAXATOMLIST); j++) - { - param.a[j] = aa[j]-1; - } - for (j = 0; (j < MAXFORCEPARAM); j++) + std::vector atoms; + for (int j = 0; (j < nral); j++) { - param.c[j] = 0.0; + atoms.emplace_back(aa[j]-1); } - + std::array forceParam = {0.0}; + InteractionType param(atoms, forceParam, ""); /* Get the cmap type for this cmap angle */ bFound = default_cmap_params(bondtype, at, atypes, ¶m, FALSE, &cmap_type, &ncmap_params, wi); @@ -2250,14 +2333,14 @@ void push_cmap(Directive d, gmx::ArrayRef bondtype, if (bFound && ncmap_params == 1) { /* Put the values in the appropriate arrays */ - param.c[0] = cmap_type; - add_param_to_list(&bond[ftype], ¶m); + param.setForceParameter(0, cmap_type); + add_param_to_list(&bond[ftype], param); } else { /* This is essentially the same check as in default_cmap_params() done one more time */ auto message = gmx::formatString("Unable to assign a cmap type to torsion %d %d %d %d and %d\n", - param.a[0]+1, param.a[1]+1, param.a[2]+1, param.a[3]+1, param.a[4]+1); + param.ai()+1, param.aj()+1, param.ak()+1, param.al()+1, param.am()+1); warning_error_and_exit(wi, message, FARGS); } } @@ -2268,22 +2351,12 @@ void push_vsitesn(Directive d, gmx::ArrayRef bond, t_atoms *at, char *line, warninp *wi) { - char *ptr; - int type, ftype, j, n, ret, nj, a; - int *atc = nullptr; - double *weight = nullptr, weight_tot; - t_param param; - - /* default force parameters */ - for (j = 0; (j < MAXATOMLIST); j++) - { - param.a[j] = NOTSET; - } - for (j = 0; (j < MAXFORCEPARAM); j++) - { - param.c[j] = 0.0; - } + char *ptr; + int type, ftype, n, ret, nj, a; + int *atc = nullptr; + double *weight = nullptr, weight_tot; + std::array forceParam = {0.0}; ptr = line; ret = sscanf(ptr, "%d%n", &a, &n); ptr += n; @@ -2294,11 +2367,10 @@ void push_vsitesn(Directive d, gmx::ArrayRef bond, warning_error_and_exit(wi, message, FARGS); } - param.a[0] = a - 1; - sscanf(ptr, "%d%n", &type, &n); ptr += n; ftype = ifunc_index(d, type); + int firstAtom = a - 1; weight_tot = 0; nj = 0; @@ -2360,13 +2432,13 @@ void push_vsitesn(Directive d, gmx::ArrayRef bond, warning_error_and_exit(wi, "The total mass of the construting atoms is zero", FARGS); } - for (j = 0; j < nj; j++) + for (int j = 0; j < nj; j++) { - param.a[1] = atc[j]; - param.c[0] = nj; - param.c[1] = weight[j]/weight_tot; + std::vector atoms = {firstAtom, atc[j]}; + forceParam[0] = nj; + forceParam[1] = weight[j]/weight_tot; /* Put the values in the appropriate arrays */ - add_param_to_list (&bond[ftype], ¶m); + add_param_to_list (&bond[ftype], InteractionType(atoms, forceParam)); } sfree(atc); @@ -2489,9 +2561,8 @@ void push_excl(char *line, gmx::ArrayRef b2, warninp *wi) int add_atomtype_decoupled(t_symtab *symtab, PreprocessingAtomTypes *at, t_nbparam ***nbparam, t_nbparam ***pair) { - t_atom atom; - t_param param; - int i, nr; + t_atom atom; + int nr; /* Define an atom type with all parameters set to zero (no interactions) */ atom.q = 0.0; @@ -2500,12 +2571,9 @@ int add_atomtype_decoupled(t_symtab *symtab, PreprocessingAtomTypes *at, * this should be changed automatically later when required. */ atom.ptype = eptAtom; - for (i = 0; (i < MAXFORCEPARAM); i++) - { - param.c[i] = 0.0; - } - nr = at->addType(symtab, atom, "decoupled", ¶m, -1, 0); + std::array forceParam = {0.0}; + nr = at->addType(symtab, atom, "decoupled", InteractionType({}, forceParam, ""), -1, 0); /* Add space in the non-bonded parameters matrix */ realloc_nb_params(at, nbparam, pair); @@ -2516,17 +2584,11 @@ int add_atomtype_decoupled(t_symtab *symtab, PreprocessingAtomTypes *at, static void convert_pairs_to_pairsQ(gmx::ArrayRef plist, real fudgeQQ, t_atoms *atoms) { - t_param *paramp1, *paramp2, *paramnew; - int i, j, p1nr, p2nr, p2newnr; - /* Add the pair list to the pairQ list */ - p1nr = plist[F_LJ14].nr; - p2nr = plist[F_LJC14_Q].nr; - p2newnr = p1nr + p2nr; - snew(paramnew, p2newnr); + std::vector paramnew; - paramp1 = plist[F_LJ14].param; - paramp2 = plist[F_LJC14_Q].param; + gmx::ArrayRef paramp1 = plist[F_LJ14].interactionTypes; + gmx::ArrayRef paramp2 = plist[F_LJC14_Q].interactionTypes; /* Fill in the new F_LJC14_Q array with the old one. NOTE: it may be possible to just ADD the converted F_LJ14 array @@ -2534,81 +2596,49 @@ static void convert_pairs_to_pairsQ(gmx::ArrayRef pli a new sized memory structure, better just to deep copy it all. */ - for (i = 0; i < p2nr; i++) - { - /* Copy over parameters */ - for (j = 0; j < 5; j++) /* entries are 0-4 for F_LJC14_Q */ - { - paramnew[i].c[j] = paramp2[i].c[j]; - } - /* copy over atoms */ - for (j = 0; j < 2; j++) - { - paramnew[i].a[j] = paramp2[i].a[j]; - } + for (const auto ¶m : paramp2) + { + paramnew.emplace_back(param); } - for (i = p2nr; i < p2newnr; i++) + for (const auto ¶m : paramp1) { - j = i-p2nr; - - /* Copy over parameters */ - paramnew[i].c[0] = fudgeQQ; - paramnew[i].c[1] = atoms->atom[paramp1[j].a[0]].q; - paramnew[i].c[2] = atoms->atom[paramp1[j].a[1]].q; - paramnew[i].c[3] = paramp1[j].c[0]; - paramnew[i].c[4] = paramp1[j].c[1]; - - /* copy over atoms */ - paramnew[i].a[0] = paramp1[j].a[0]; - paramnew[i].a[1] = paramp1[j].a[1]; + std::vector forceParam = { + fudgeQQ, atoms->atom[param.ai()].q, atoms->atom[param.aj()].q, + param.c0(), param.c1() + }; + paramnew.emplace_back(InteractionType(param.atoms(), forceParam, "")); } - /* free the old pairlists */ - sfree(plist[F_LJC14_Q].param); - sfree(plist[F_LJ14].param); - /* now assign the new data to the F_LJC14_Q structure */ - plist[F_LJC14_Q].param = paramnew; - plist[F_LJC14_Q].nr = p2newnr; + plist[F_LJC14_Q].interactionTypes = paramnew; /* Empty the LJ14 pairlist */ - plist[F_LJ14].nr = 0; - plist[F_LJ14].param = nullptr; + plist[F_LJ14].interactionTypes.clear(); } static void generate_LJCpairsNB(MoleculeInformation *mol, int nb_funct, InteractionTypeParameters *nbp, warninp *wi) { - int n, ntype, i, j, k; - t_atom *atom; - t_blocka *excl; - bool bExcl; - t_param param; + int n, ntype; + t_atom *atom; + t_blocka *excl; + bool bExcl; n = mol->atoms.nr; atom = mol->atoms.atom; - ntype = static_cast(std::sqrt(static_cast(nbp->nr))); - GMX_ASSERT(ntype * ntype == nbp->nr, "Number of pairs of generated non-bonded parameters should be a perfect square"); - - for (i = 0; i < MAXATOMLIST; i++) - { - param.a[i] = NOTSET; - } - for (i = 0; i < MAXFORCEPARAM; i++) - { - param.c[i] = NOTSET; - } + ntype = static_cast(std::sqrt(static_cast(nbp->size()))); + GMX_ASSERT(ntype * ntype == gmx::ssize(*nbp), "Number of pairs of generated non-bonded parameters should be a perfect square"); /* Add a pair interaction for all non-excluded atom pairs */ excl = &mol->excls; - for (i = 0; i < n; i++) + for (int i = 0; i < n; i++) { - for (j = i+1; j < n; j++) + for (int j = i+1; j < n; j++) { bExcl = FALSE; - for (k = excl->index[i]; k < excl->index[i+1]; k++) + for (int k = excl->index[i]; k < excl->index[i+1]; k++) { if (excl->a[k] == j) { @@ -2624,13 +2654,14 @@ static void generate_LJCpairsNB(MoleculeInformation *mol, int nb_funct, Interact "for Van der Waals type Lennard-Jones"); warning_error_and_exit(wi, message, FARGS); } - param.a[0] = i; - param.a[1] = j; - param.c[0] = atom[i].q; - param.c[1] = atom[j].q; - param.c[2] = nbp->param[ntype*atom[i].type+atom[j].type].c[0]; - param.c[3] = nbp->param[ntype*atom[i].type+atom[j].type].c[1]; - add_param_to_list(&mol->plist[F_LJC_PAIRS_NB], ¶m); + std::vector atoms = {i, j}; + std::vector forceParam = { + atom[i].q, + atom[j].q, + nbp->interactionTypes[ntype*atom[i].type+atom[j].type].c0(), + nbp->interactionTypes[ntype*atom[i].type+atom[j].type].c1() + }; + add_param_to_list(&mol->plist[F_LJC_PAIRS_NB], InteractionType(atoms, forceParam)); } } } diff --git a/src/gromacs/gmxpreprocess/toppush.h b/src/gromacs/gmxpreprocess/toppush.h index 33371cafc7..aa0a3441f2 100644 --- a/src/gromacs/gmxpreprocess/toppush.h +++ b/src/gromacs/gmxpreprocess/toppush.h @@ -50,7 +50,7 @@ struct t_atoms; struct t_block; struct MoleculeInformation; struct t_nbparam; -struct t_param; +class InteractionType; struct InteractionTypeParameters; struct PreprocessResidue; struct warninp; diff --git a/src/gromacs/gmxpreprocess/topshake.cpp b/src/gromacs/gmxpreprocess/topshake.cpp index ca0b56ff7b..5129535d42 100644 --- a/src/gromacs/gmxpreprocess/topshake.cpp +++ b/src/gromacs/gmxpreprocess/topshake.cpp @@ -53,35 +53,7 @@ #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/smalloc.h" -static void copy_bond (InteractionTypeParameters *pr, int to, int from) -/* copies an entry in a bond list to another position. - * does no allocing or freeing of memory - */ -{ - /*memcpy((char*) &(pr->param[to]),(char*) &(pr->param[from]), - (size_t)sizeof(pr->param[from]));*/ - int i; - - if (to != from) - { - range_check(to, 0, pr->nr); - range_check(from, 0, pr->nr); - for (i = 0; (i < MAXATOMLIST); i++) - { - pr->param[to].a[i] = pr->param[from].a[i]; - } - for (i = 0; (i < MAXFORCEPARAM); i++) - { - pr->param[to].c[i] = pr->param[from].c[i]; - } - for (i = 0; (i < MAXSLEN); i++) - { - pr->param[to].s[i] = pr->param[from].s[i]; - } - } -} - -static int count_hydrogens (char ***atomname, int nra, const int a[]) +static int count_hydrogens (char ***atomname, int nra, gmx::ArrayRef a) { int nh; @@ -107,8 +79,6 @@ void make_shake(gmx::ArrayRef plist, t_atoms *atoms, { char ***info = atoms->atomname; real b_ij, b_jk; - int i, j; - if (nshake != eshNONE) { switch (nshake) @@ -140,72 +110,69 @@ void make_shake(gmx::ArrayRef plist, t_atoms *atoms, { InteractionTypeParameters *bonds = &(plist[ftype]); - for (int ftype_a = 0; (bonds->nr > 0 && ftype_a < F_NRE); ftype_a++) + for (int ftype_a = 0; (gmx::ssize(*bonds) > 0 && ftype_a < F_NRE); ftype_a++) { if (interaction_function[ftype_a].flags & IF_ATYPE) { InteractionTypeParameters *pr = &(plist[ftype_a]); - for (int i = 0; (i < pr->nr); ) + for (auto parm = pr->interactionTypes.begin(); parm != pr->interactionTypes.end(); ) { - t_param *ang = &(pr->param[i]); + const InteractionType *ang = &(*parm); #ifdef DEBUG printf("Angle: %d-%d-%d\n", ang->ai(), ang->aj(), ang->ak()); #endif - int numhydrogens = count_hydrogens(info, 3, ang->a); + int numhydrogens = count_hydrogens(info, 3, ang->atoms()); if ((nshake == eshALLANGLES) || (numhydrogens > 1) || - (numhydrogens == 1 && toupper(**(info[ang->a[1]])) == 'O')) + (numhydrogens == 1 && toupper(**(info[ang->aj()])) == 'O')) { /* Can only add hydrogen angle shake, if the two bonds * are constrained. * append this angle to the shake list */ - t_param p; - p.ai() = ang->ai(); - p.aj() = ang->ak(); + std::vector atomNumbers = {ang->ai(), ang->ak()}; /* Calculate length of constraint */ bool bFound = false; b_ij = b_jk = 0.0; - for (int j = 0; !bFound && (j < bonds->nr); j++) + for (const auto &bond : bonds->interactionTypes) { - t_param *bond = &(bonds->param[j]); - if (((bond->ai() == ang->ai()) && - (bond->aj() == ang->aj())) || - ((bond->ai() == ang->aj()) && - (bond->aj() == ang->ai()))) + if (((bond.ai() == ang->ai()) && + (bond.aj() == ang->aj())) || + ((bond.ai() == ang->aj()) && + (bond.aj() == ang->ai()))) { - b_ij = bond->c0(); + b_ij = bond.c0(); } - if (((bond->ai() == ang->ak()) && - (bond->aj() == ang->aj())) || - ((bond->ai() == ang->aj()) && - (bond->aj() == ang->ak()))) + if (((bond.ai() == ang->ak()) && + (bond.aj() == ang->aj())) || + ((bond.ai() == ang->aj()) && + (bond.aj() == ang->ak()))) { - b_jk = bond->c0(); + b_jk = bond.c0(); } bFound = (b_ij != 0.0) && (b_jk != 0.0); } if (bFound) { + real param = std::sqrt( b_ij*b_ij + b_jk*b_jk - + 2.0*b_ij*b_jk*cos(DEG2RAD*ang->c0())); + std::vector forceParm = {param, param}; /* apply law of cosines */ - p.c0() = std::sqrt( b_ij*b_ij + b_jk*b_jk - - 2.0*b_ij*b_jk*cos(DEG2RAD*ang->c0()) ); - p.c1() = p.c0(); #ifdef DEBUG - printf("p: %d, q: %d, dist: %12.5e\n", p.ai(), p.aj(), p.c0()); + printf("p: %d, q: %d, dist: %12.5e\n", atomNumbers[0], + atomNumbers[1], forceParm[0]); #endif - add_param_to_list (&(plist[F_CONSTR]), &p); + add_param_to_list (&(plist[F_CONSTR]), InteractionType(atomNumbers, forceParm)); /* move the last bond to this position */ - copy_bond (pr, i, pr->nr-1); - /* should free memory here!! */ - pr->nr--; + *parm = *(pr->interactionTypes.end() - 1); + pr->interactionTypes.erase(pr->interactionTypes.end() - 1); } } else { - i++; + ++parm; } } } /* if IF_ATYPE */ @@ -222,26 +189,22 @@ void make_shake(gmx::ArrayRef plist, t_atoms *atoms, if (interaction_function[ftype].flags & IF_BTYPE) { InteractionTypeParameters *pr = &(plist[ftype]); - j = 0; - for (i = 0; i < pr->nr; i++) + for (auto parm = pr->interactionTypes.begin(); parm != pr->interactionTypes.end(); ) { if ( (nshake != eshHBONDS) || - (count_hydrogens (info, 2, pr->param[i].a) > 0) ) + (count_hydrogens (info, 2, parm->atoms()) > 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(); - p.c1() = pr->param[i].c2(); - add_param_to_list (&(plist[F_CONSTR]), &p); + std::vector atomNumbers = {parm->ai(), parm->aj()}; + std::vector forceParm = { parm->c0(), parm->c2()}; + add_param_to_list (&(plist[F_CONSTR]), InteractionType(atomNumbers, forceParm)); + parm = pr->interactionTypes.erase(parm); } else { - copy_bond(pr, j++, i); + ++parm; } } - pr->nr = j; } } } diff --git a/src/gromacs/gmxpreprocess/toputil.cpp b/src/gromacs/gmxpreprocess/toputil.cpp index 5fc3e04b0f..8a1dd69412 100644 --- a/src/gromacs/gmxpreprocess/toputil.cpp +++ b/src/gromacs/gmxpreprocess/toputil.cpp @@ -57,112 +57,9 @@ /* UTILITIES */ -void set_p_string(t_param *p, const char *s) +void add_param_to_list(InteractionTypeParameters *list, const InteractionType &b) { - if (s) - { - if (strlen(s) < sizeof(p->s)-1) - { - strncpy(p->s, s, sizeof(p->s)); - } - else - { - gmx_fatal(FARGS, "Increase MAXSLEN in the grompp code to at least %zu," - " or shorten your definition of bonds like %s to at most %d", - strlen(s)+1, s, MAXSLEN-1); - } - } - else - { - strcpy(p->s, ""); - } -} - -void pr_alloc (int extra, InteractionTypeParameters *pr) -{ - int i, j; - - /* get new space for arrays */ - if (extra < 0) - { - gmx_fatal(FARGS, "Trying to make array smaller.\n"); - } - if (extra == 0) - { - return; - } - GMX_ASSERT(pr->nr != 0 || pr->param == nullptr, "Invalid InteractionTypeParameters object"); - if (pr->nr+extra > pr->maxnr) - { - pr->maxnr = std::max(static_cast(1.2*pr->maxnr), pr->maxnr + extra); - srenew(pr->param, pr->maxnr); - for (i = pr->nr; (i < pr->maxnr); i++) - { - for (j = 0; (j < MAXATOMLIST); j++) - { - pr->param[i].a[j] = 0; - } - for (j = 0; (j < MAXFORCEPARAM); j++) - { - pr->param[i].c[j] = 0; - } - set_p_string(&(pr->param[i]), ""); - } - } -} - -void init_plist(gmx::ArrayRef plist) -{ - int i; - - for (i = 0; (i < F_NRE); i++) - { - plist[i].nr = 0; - plist[i].maxnr = 0; - plist[i].param = nullptr; - } -} - -void done_plist(gmx::ArrayRef plist) -{ - for (int i = 0; i < F_NRE; i++) - { - sfree(plist[i].param); - } -} - -void cp_param(t_param *dest, const t_param *src) -{ - int j; - - for (j = 0; (j < MAXATOMLIST); j++) - { - dest->a[j] = src->a[j]; - } - for (j = 0; (j < MAXFORCEPARAM); j++) - { - dest->c[j] = src->c[j]; - } - strncpy(dest->s, src->s, sizeof(dest->s)); -} - -void add_param_to_list(InteractionTypeParameters *list, t_param *b) -{ - /* allocate one position extra */ - pr_alloc (1, list); - - /* fill the arrays */ - for (int j = 0; (j < MAXFORCEPARAM); j++) - { - list->param[list->nr].c[j] = b->c[j]; - } - for (int j = 0; (j < MAXATOMLIST); j++) - { - list->param[list->nr].a[j] = b->a[j]; - } - memset(list->param[list->nr].s, 0, sizeof(list->param[list->nr].s)); - - list->nr++; + list->interactionTypes.emplace_back(b); } /* PRINTING STRUCTURES */ @@ -180,7 +77,7 @@ static void print_bt(FILE *out, Directive d, PreprocessingAtomTypes *at, const InteractionTypeParameters *bt = &(plist[ftype]); - if (bt->nr == 0) + if (bt->size() == 0) { return; } @@ -285,34 +182,36 @@ static void print_bt(FILE *out, Directive d, PreprocessingAtomTypes *at, fprintf (out, "\n"); /* print bondtypes */ - for (int i = 0; (i < bt->nr); i++) + for (const auto &parm : bt->interactionTypes) { - bSwapParity = (bt->param[i].c0() == NOTSET) && (bt->param[i].c1() == -1); + bSwapParity = (parm.c0() == NOTSET) && (parm.c1() == -1); + gmx::ArrayRef atoms = parm.atoms(); if (!bDih) { for (int j = 0; (j < nral); j++) { - fprintf (out, "%5s ", at->atomNameFromAtomType(bt->param[i].a[j])); + fprintf (out, "%5s ", at->atomNameFromAtomType(atoms[j])); } } else { for (int j = 0; (j < 2); j++) { - fprintf (out, "%5s ", at->atomNameFromAtomType(bt->param[i].a[dihp[f][j]])); + fprintf (out, "%5s ", at->atomNameFromAtomType(atoms[dihp[f][j]])); } } fprintf (out, "%5d ", bSwapParity ? -f-1 : f+1); - if (bt->param[i].s[0]) + if (!parm.interactionTypeName().empty()) { - fprintf(out, " %s", bt->param[i].s); + fprintf(out, " %s", parm.interactionTypeName().c_str()); } else { - for (int j = 0; (j < nrfp && (bt->param[i].c[j] != NOTSET)); j++) + gmx::ArrayRef forceParam = parm.forceParam(); + for (int j = 0; (j < nrfp) && (forceParam[j] != NOTSET); j++) { - fprintf (out, "%13.6e ", bt->param[i].c[j]); + fprintf (out, "%13.6e ", forceParam[j]); } } @@ -490,22 +389,19 @@ void print_bondeds(FILE *out, int natoms, Directive d, int ftype, int fsubtype, gmx::ArrayRef plist) { t_symtab stab; - t_param *param; t_atom *a; PreprocessingAtomTypes atype; snew(a, 1); - snew(param, 1); open_symtab(&stab); for (int i = 0; (i < natoms); i++) { char buf[12]; sprintf(buf, "%4d", (i+1)); - atype.addType(&stab, *a, buf, param, 0, 0); + atype.addType(&stab, *a, buf, InteractionType({}, {}), 0, 0); } print_bt(out, d, &atype, ftype, fsubtype, plist, TRUE); done_symtab(&stab); sfree(a); - sfree(param); } diff --git a/src/gromacs/gmxpreprocess/toputil.h b/src/gromacs/gmxpreprocess/toputil.h index 2891177b17..1c04e64704 100644 --- a/src/gromacs/gmxpreprocess/toputil.h +++ b/src/gromacs/gmxpreprocess/toputil.h @@ -50,25 +50,15 @@ struct t_atoms; struct t_blocka; struct t_excls; struct MoleculeInformation; -struct t_param; +class InteractionType; struct InteractionTypeParameters; /* UTILITIES */ int name2index(char *str, char ***typenames, int ntypes); -void pr_alloc (int extra, InteractionTypeParameters *pr); +void add_param_to_list(InteractionTypeParameters *list, const InteractionType &b); -void set_p_string(t_param *p, const char *s); - -void cp_param(t_param *dest, const t_param *src); - -void add_param_to_list(InteractionTypeParameters *list, t_param *b); - -/* INITIATE */ - -void init_plist(gmx::ArrayRef plist); -void done_plist(gmx::ArrayRef plist); /* PRINTING */ diff --git a/src/gromacs/gmxpreprocess/vsite_parm.cpp b/src/gromacs/gmxpreprocess/vsite_parm.cpp index c0c3b36dbd..88635cc9c1 100644 --- a/src/gromacs/gmxpreprocess/vsite_parm.cpp +++ b/src/gromacs/gmxpreprocess/vsite_parm.cpp @@ -73,17 +73,21 @@ typedef struct { t_iatom &al() { return a[3]; } } t_mybonded; -typedef struct { - int ftype; - t_param *param; -} vsitebondparam_t; +struct VsiteBondParameter +{ + VsiteBondParameter(int ftype, const InteractionType &type) + : ftype_(ftype), type_(type) + {} + int ftype_; + const InteractionType &type_; +}; struct Atom2VsiteBond { //! Function type for conversion. - int ftype; + int ftype; //! The vsite parameters in a list. - std::vector vSiteBondedParameters; + std::vector vSiteBondedParameters; }; typedef struct { @@ -108,24 +112,24 @@ static int vsite_bond_nrcheck(int ftype) } static void enter_bonded(int nratoms, int *nrbonded, t_mybonded **bondeds, - const t_param *param) + const InteractionType &type) { - int j; - srenew(*bondeds, *nrbonded+1); /* copy atom numbers */ - for (j = 0; j < nratoms; j++) + gmx::ArrayRef atoms = type.atoms(); + GMX_RELEASE_ASSERT(nratoms == atoms.ssize(), "Size of atom array must much"); + for (int j = 0; j < nratoms; j++) { - (*bondeds)[*nrbonded].a[j] = param->a[j]; + (*bondeds)[*nrbonded].a[j] = atoms[j]; } /* copy parameter */ - (*bondeds)[*nrbonded].c = param->c0(); + (*bondeds)[*nrbonded].c = type.c0(); (*nrbonded)++; } -static void get_bondeds(int nrat, const t_iatom atoms[], +static void get_bondeds(int nrat, gmx::ArrayRef atoms, gmx::ArrayRef at2vb, int *nrbond, t_mybonded **bonds, int *nrang, t_mybonded **angles, @@ -135,15 +139,15 @@ static void get_bondeds(int nrat, const t_iatom atoms[], { for (auto &vsite : at2vb[atoms[k]].vSiteBondedParameters) { - int ftype = vsite.ftype; - const t_param *param = vsite.param; - int nrcheck = vsite_bond_nrcheck(ftype); + int ftype = vsite.ftype_; + const InteractionType &type = vsite.type_; + int nrcheck = vsite_bond_nrcheck(ftype); /* abuse nrcheck to see if we're adding bond, angle or idih */ switch (nrcheck) { - case 2: enter_bonded(nrcheck, nrbond, bonds, param); break; - case 3: enter_bonded(nrcheck, nrang, angles, param); break; - case 4: enter_bonded(nrcheck, nridih, idihs, param); break; + case 2: enter_bonded(nrcheck, nrbond, bonds, type); break; + case 3: enter_bonded(nrcheck, nrang, angles, type); break; + case 4: enter_bonded(nrcheck, nridih, idihs, type); break; } } } @@ -153,7 +157,6 @@ static std::vector make_at2vsitebond(int natoms, gmx::ArrayRef plist) { bool *bVSI; - t_iatom *aa; std::vector at2vb(natoms); @@ -162,11 +165,12 @@ make_at2vsitebond(int natoms, gmx::ArrayRef plist) { if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN) { - for (int i = 0; (i < plist[ftype].nr); i++) + for (int i = 0; (i < gmx::ssize(plist[ftype])); i++) { + gmx::ArrayRef atoms = plist[ftype].interactionTypes[i].atoms(); for (int j = 0; j < NRAL(ftype); j++) { - bVSI[plist[ftype].param[i].a[j]] = TRUE; + bVSI[atoms[j]] = TRUE; } } } @@ -177,16 +181,14 @@ make_at2vsitebond(int natoms, gmx::ArrayRef plist) int nrcheck = vsite_bond_nrcheck(ftype); if (nrcheck > 0) { - for (int i = 0; (i < plist[ftype].nr); i++) + for (int i = 0; (i < gmx::ssize(plist[ftype])); i++) { - aa = plist[ftype].param[i].a; + gmx::ArrayRef aa = plist[ftype].interactionTypes[i].atoms(); for (int j = 0; j < nrcheck; j++) { if (bVSI[aa[j]]) { - at2vb[aa[j]].vSiteBondedParameters.emplace_back(); - at2vb[aa[j]].vSiteBondedParameters.back().ftype = ftype; - at2vb[aa[j]].vSiteBondedParameters.back().param = &plist[ftype].param[i]; + at2vb[aa[j]].vSiteBondedParameters.emplace_back(ftype, plist[ftype].interactionTypes[i]); } } } @@ -200,38 +202,38 @@ make_at2vsitebond(int natoms, gmx::ArrayRef plist) static at2vsitecon_t *make_at2vsitecon(int natoms, gmx::ArrayRef plist) { bool *bVSI; - int ftype, i, j, ai, aj, nr; at2vsitecon_t *at2vc; snew(at2vc, natoms); snew(bVSI, natoms); - for (ftype = 0; (ftype < F_NRE); ftype++) + for (int ftype = 0; (ftype < F_NRE); ftype++) { if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN) { - for (i = 0; (i < plist[ftype].nr); i++) + for (int i = 0; (i < gmx::ssize(plist[ftype])); i++) { - for (j = 0; j < NRAL(ftype); j++) + gmx::ArrayRef atoms = plist[ftype].interactionTypes[i].atoms(); + for (int j = 0; j < NRAL(ftype); j++) { - bVSI[plist[ftype].param[i].a[j]] = TRUE; + bVSI[atoms[j]] = TRUE; } } } } - for (ftype = 0; (ftype < F_NRE); ftype++) + for (int ftype = 0; (ftype < F_NRE); ftype++) { if (interaction_function[ftype].flags & IF_CONSTRAINT) { - for (i = 0; (i < plist[ftype].nr); i++) + for (int i = 0; (i < gmx::ssize(plist[ftype])); i++) { - ai = plist[ftype].param[i].ai(); - aj = plist[ftype].param[i].aj(); + int ai = plist[ftype].interactionTypes[i].ai(); + int aj = plist[ftype].interactionTypes[i].aj(); if (bVSI[ai] && bVSI[aj]) { /* Store forward direction */ - nr = at2vc[ai].nr; + int nr = at2vc[ai].nr; if (nr % 10 == 0) { srenew(at2vc[ai].aj, nr+10); @@ -257,9 +259,7 @@ static at2vsitecon_t *make_at2vsitecon(int natoms, gmx::ArrayRef forceParam = type.forceParam(); + for (int j = 0; j < NRFP(ftype); j++) { - fprintf(fp, ".c[%d]=%g ", j, param->c[j]); + fprintf(fp, ".c[%d]=%g ", j, forceParam[j]); } fprintf(fp, "\n"); pass++; @@ -395,7 +393,7 @@ static const char *get_atomtype_name_AB(t_atom *atom, PreprocessingAtomTypes *at } static bool calc_vsite3_param(PreprocessingAtomTypes *atypes, - t_param *param, t_atoms *at, + InteractionType *param, t_atoms *at, int nrbond, t_mybonded *bonds, int nrang, t_mybonded *angles ) { @@ -467,14 +465,13 @@ static bool calc_vsite3_param(PreprocessingAtomTypes *atypes, gmx_fatal(FARGS, "calc_vsite3_param not implemented for the general case " "(atom %d)", param->ai()+1); } - - param->c0() = a; - param->c1() = b; + param->setForceParameter(0, a); + param->setForceParameter(1, b); return bError; } -static bool calc_vsite3fd_param(t_param *param, +static bool calc_vsite3fd_param(InteractionType *param, int nrbond, t_mybonded *bonds, int nrang, t_mybonded *angles) { @@ -496,13 +493,13 @@ static bool calc_vsite3fd_param(t_param *param, rk = bjk * std::sin(aijk); rl = bjl * std::sin(aijl); - param->c0() = rk / (rk + rl); - param->c1() = -bij; /* 'bond'-length for fixed distance vsite */ + param->setForceParameter(0, rk / (rk + rl)); + param->setForceParameter(1, -bij); return bError; } -static bool calc_vsite3fad_param(t_param *param, +static bool calc_vsite3fad_param(InteractionType *param, int nrbond, t_mybonded *bonds, int nrang, t_mybonded *angles) { @@ -521,19 +518,19 @@ static bool calc_vsite3fad_param(t_param *param, aijk = get_angle (nrang, angles, param->ai(), param->aj(), param->ak()); bError = (bij == NOTSET) || (aijk == NOTSET); - param->c1() = bij; /* 'bond'-length for fixed distance vsite */ - param->c0() = RAD2DEG*aijk; /* 'bond'-angle for fixed angle vsite */ + param->setForceParameter(1, bij); + param->setForceParameter(0, RAD2DEG*aijk); if (bSwapParity) { - param->c0() = 360 - param->c0(); + param->setForceParameter(0, 360 - param->c0()); } return bError; } static bool calc_vsite3out_param(PreprocessingAtomTypes *atypes, - t_param *param, t_atoms *at, + InteractionType *param, t_atoms *at, int nrbond, t_mybonded *bonds, int nrang, t_mybonded *angles) { @@ -616,20 +613,20 @@ static bool calc_vsite3out_param(PreprocessingAtomTypes *atypes, / ( bjk*bjl*std::sin(akjl) ); } - param->c0() = a; - param->c1() = b; + param->setForceParameter(0, a); + param->setForceParameter(1, b); if (bSwapParity) { - param->c2() = -c; + param->setForceParameter(2, -c); } else { - param->c2() = c; + param->setForceParameter(2, c); } return bError; } -static bool calc_vsite4fd_param(t_param *param, +static bool calc_vsite4fd_param(InteractionType *param, int nrbond, t_mybonded *bonds, int nrang, t_mybonded *angles) { @@ -676,9 +673,9 @@ static bool calc_vsite4fd_param(t_param *param, cl = -pk / ( pl*cosakl - pk + pl*sinakl*(pm*cosakm-pk)/(pm*sinakm) ); cm = -pk / ( pm*cosakm - pk + pm*sinakm*(pl*cosakl-pk)/(pl*sinakl) ); - param->c0() = cl; - param->c1() = cm; - param->c2() = -bij; + param->setForceParameter(0, cl); + param->setForceParameter(1, cm); + param->setForceParameter(2, -bij); } return bError; @@ -686,7 +683,7 @@ static bool calc_vsite4fd_param(t_param *param, static bool -calc_vsite4fdn_param(t_param *param, +calc_vsite4fdn_param(InteractionType *param, int nrbond, t_mybonded *bonds, int nrang, t_mybonded *angles) { @@ -733,10 +730,9 @@ calc_vsite4fdn_param(t_param *param, a = pk/pl; b = pk/pm; - param->c0() = a; - param->c1() = b; - param->c2() = bij; - + param->setForceParameter(0, a); + param->setForceParameter(1, b); + param->setForceParameter(2, bij); } return bError; @@ -747,9 +743,9 @@ calc_vsite4fdn_param(t_param *param, int set_vsites(bool bVerbose, t_atoms *atoms, PreprocessingAtomTypes *atypes, gmx::ArrayRef plist) { - int i, j, ftype; + int ftype; int nvsite, nrbond, nrang, nridih, nrset; - bool bFirst, bSet, bERROR; + bool bFirst, bERROR; t_mybonded *bonds; t_mybonded *angles; t_mybonded *idihs; @@ -764,7 +760,7 @@ int set_vsites(bool bVerbose, t_atoms *atoms, PreprocessingAtomTypes *atypes, { if (interaction_function[ftype].flags & IF_VSITE) { - nvsite += plist[ftype].nr; + nvsite += plist[ftype].size(); if (ftype == F_VSITEN) { @@ -773,19 +769,21 @@ int set_vsites(bool bVerbose, t_atoms *atoms, PreprocessingAtomTypes *atypes, } nrset = 0; - for (i = 0; (i < plist[ftype].nr); i++) + int i = 0; + for (auto ¶m : plist[ftype].interactionTypes) { /* check if all parameters are set */ - bSet = TRUE; - for (j = 0; j < NRFP(ftype) && bSet; j++) + bool bSet = true; + gmx::ArrayRef forceParam = param.forceParam(); + for (int j = 0; (j < NRFP(ftype)) && bSet; j++) { - bSet = plist[ftype].param[i].c[j] != NOTSET; + bSet = forceParam[j] != NOTSET; } if (debug) { fprintf(debug, "bSet=%s ", gmx::boolToString(bSet)); - print_param(debug, ftype, i, &plist[ftype].param[i]); + printInteractionType(debug, ftype, i, plist[ftype].interactionTypes[i]); } if (!bSet) { @@ -801,13 +799,13 @@ int set_vsites(bool bVerbose, t_atoms *atoms, PreprocessingAtomTypes *atypes, idihs = nullptr; nrset++; /* now set the vsite parameters: */ - get_bondeds(NRAL(ftype), plist[ftype].param[i].a, at2vb, + get_bondeds(NRAL(ftype), param.atoms(), at2vb, &nrbond, &bonds, &nrang, &angles, &nridih, &idihs); if (debug) { fprintf(debug, "Found %d bonds, %d angles and %d idihs " "for virtual site %d (%s)\n", nrbond, nrang, nridih, - plist[ftype].param[i].ai()+1, + param.ai()+1, interaction_function[ftype].longname); print_bad(debug, nrbond, bonds, nrang, angles, nridih, idihs); } /* debug */ @@ -815,39 +813,39 @@ int set_vsites(bool bVerbose, t_atoms *atoms, PreprocessingAtomTypes *atypes, { case F_VSITE3: bERROR = - calc_vsite3_param(atypes, &(plist[ftype].param[i]), atoms, + calc_vsite3_param(atypes, ¶m, atoms, nrbond, bonds, nrang, angles); break; case F_VSITE3FD: bERROR = - calc_vsite3fd_param(&(plist[ftype].param[i]), + calc_vsite3fd_param(¶m, nrbond, bonds, nrang, angles); break; case F_VSITE3FAD: bERROR = - calc_vsite3fad_param(&(plist[ftype].param[i]), + calc_vsite3fad_param(¶m, nrbond, bonds, nrang, angles); break; case F_VSITE3OUT: bERROR = - calc_vsite3out_param(atypes, &(plist[ftype].param[i]), atoms, + calc_vsite3out_param(atypes, ¶m, atoms, nrbond, bonds, nrang, angles); break; case F_VSITE4FD: bERROR = - calc_vsite4fd_param(&(plist[ftype].param[i]), + calc_vsite4fd_param(¶m, nrbond, bonds, nrang, angles); break; case F_VSITE4FDN: bERROR = - calc_vsite4fdn_param(&(plist[ftype].param[i]), + calc_vsite4fdn_param(¶m, nrbond, bonds, nrang, angles); break; default: gmx_fatal(FARGS, "Automatic parameter generation not supported " "for %s atom %d", interaction_function[ftype].longname, - plist[ftype].param[i].ai()+1); + param.ai()+1); bERROR = TRUE; } /* switch */ if (bERROR) @@ -855,13 +853,14 @@ int set_vsites(bool bVerbose, t_atoms *atoms, PreprocessingAtomTypes *atypes, gmx_fatal(FARGS, "Automatic parameter generation not supported " "for %s atom %d for this bonding configuration", interaction_function[ftype].longname, - plist[ftype].param[i].ai()+1); + param.ai()+1); } sfree(bonds); sfree(angles); sfree(idihs); } /* if bSet */ - } /* for i */ + i++; + } } /* if IF_VSITE */ } @@ -911,21 +910,17 @@ typedef struct { static void check_vsite_constraints(gmx::ArrayRef plist, int cftype, const int vsite_type[]) { - int i, k, n; - int atom; - InteractionTypeParameters *ps; - - n = 0; - ps = &(plist[cftype]); - for (i = 0; (i < ps->nr); i++) + int n = 0; + for (const auto ¶m : plist[cftype].interactionTypes) { - for (k = 0; k < 2; k++) + gmx::ArrayRef atoms = param.atoms(); + for (int k = 0; k < 2; k++) { - atom = ps->param[i].a[k]; + int atom = atoms[k]; if (vsite_type[atom] != NOTSET) { fprintf(stderr, "ERROR: Cannot have constraint (%d-%d) with virtual site (%d)\n", - ps->param[i].ai()+1, ps->param[i].aj()+1, atom+1); + param.ai()+1, param.aj()+1, atom+1); n++; } } @@ -939,10 +934,10 @@ static void check_vsite_constraints(gmx::ArrayRef pli static void clean_vsite_bonds(gmx::ArrayRef plist, t_pindex pindex[], int cftype, const int vsite_type[]) { - int ftype, i, j, k, m, n, nvsite, nOut, kept_i; + int ftype, nOut; int nconverted, nremoved; - int atom, oatom, at1, at2; - bool bKeep, bRemove, bUsed, bPresent, bThisFD, bThisOUT, bAllFD, bFirstTwo; + int oatom, at1, at2; + bool bKeep, bRemove, bAllFD; InteractionTypeParameters *ps; if (cftype == F_CONNBONDS) @@ -951,42 +946,42 @@ static void clean_vsite_bonds(gmx::ArrayRef plist, t_ } ps = &(plist[cftype]); - kept_i = 0; nconverted = 0; nremoved = 0; nOut = 0; - for (i = 0; (i < ps->nr); i++) /* for all bonds in the plist */ + for (auto bond = ps->interactionTypes.begin(); bond != ps->interactionTypes.end(); ) { int vsnral = 0; const int *first_atoms = nullptr; - bKeep = FALSE; - bRemove = FALSE; - bAllFD = TRUE; + bKeep = false; + bRemove = false; + bAllFD = true; /* check if all virtual sites are constructed from the same atoms */ - nvsite = 0; - for (k = 0; (k < 2) && !bKeep && !bRemove; k++) + int nvsite = 0; + gmx::ArrayRef atoms = bond->atoms(); + for (int k = 0; (k < 2) && !bKeep && !bRemove; k++) { /* for all atoms in the bond */ - atom = ps->param[i].a[k]; + int atom = atoms[k]; if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN) { nvsite++; - bThisFD = ( (pindex[atom].ftype == F_VSITE3FD ) || - (pindex[atom].ftype == F_VSITE3FAD) || - (pindex[atom].ftype == F_VSITE4FD ) || - (pindex[atom].ftype == F_VSITE4FDN ) ); - bThisOUT = ( (pindex[atom].ftype == F_VSITE3OUT) && - ((interaction_function[cftype].flags & IF_CONSTRAINT) != 0u) ); + bool bThisFD = ( (pindex[atom].ftype == F_VSITE3FD ) || + (pindex[atom].ftype == F_VSITE3FAD) || + (pindex[atom].ftype == F_VSITE4FD ) || + (pindex[atom].ftype == F_VSITE4FDN ) ); + bool bThisOUT = ( (pindex[atom].ftype == F_VSITE3OUT) && + ((interaction_function[cftype].flags & IF_CONSTRAINT) != 0u) ); bAllFD = bAllFD && bThisFD; if (bThisFD || bThisOUT) { - oatom = ps->param[i].a[1-k]; /* the other atom */ + oatom = atoms[1-k]; /* the other atom */ if (vsite_type[oatom] == NOTSET && - oatom == plist[pindex[atom].ftype].param[pindex[atom].parnr].aj()) + oatom == plist[pindex[atom].ftype].interactionTypes[pindex[atom].parnr].aj()) { /* if the other atom isn't a vsite, and it is AI */ - bRemove = TRUE; + bRemove = true; if (bThisOUT) { nOut++; @@ -1006,7 +1001,7 @@ static void clean_vsite_bonds(gmx::ArrayRef plist, t_ a C++ "vector view" class" with an STL-container-like interface. */ vsnral = NRAL(pindex[atom].ftype) - 1; - first_atoms = plist[pindex[atom].ftype].param[pindex[atom].parnr].a + 1; + first_atoms = plist[pindex[atom].ftype].interactionTypes[pindex[atom].parnr].atoms().data() + 1; } else { @@ -1016,28 +1011,28 @@ static void clean_vsite_bonds(gmx::ArrayRef plist, t_ check if this vsite is constructed from the same atoms */ if (vsnral == NRAL(pindex[atom].ftype)-1) { - for (m = 0; (m < vsnral) && !bKeep; m++) + for (int m = 0; (m < vsnral) && !bKeep; m++) { const int *atoms; - bPresent = FALSE; - atoms = plist[pindex[atom].ftype].param[pindex[atom].parnr].a + 1; - for (n = 0; (n < vsnral) && !bPresent; n++) + bool bPresent = false; + atoms = plist[pindex[atom].ftype].interactionTypes[pindex[atom].parnr].atoms().data() + 1; + for (int n = 0; (n < vsnral) && !bPresent; n++) { if (atoms[m] == first_atoms[n]) { - bPresent = TRUE; + bPresent = true; } } if (!bPresent) { - bKeep = TRUE; + bKeep = true; } } } else { - bKeep = TRUE; + bKeep = true; } } } @@ -1046,40 +1041,40 @@ static void clean_vsite_bonds(gmx::ArrayRef plist, t_ if (bRemove) { - bKeep = FALSE; + bKeep = false; } else { /* if we have no virtual sites in this bond, keep it */ if (nvsite == 0) { - bKeep = TRUE; + bKeep = true; } /* TODO This loop and the corresponding loop in check_vsite_angles should be refactored into a common function */ /* check if all non-vsite atoms are used in construction: */ - bFirstTwo = TRUE; - for (k = 0; (k < 2) && !bKeep; k++) /* for all atoms in the bond */ + bool bFirstTwo = true; + for (int k = 0; (k < 2) && !bKeep; k++) /* for all atoms in the bond */ { - atom = ps->param[i].a[k]; + int atom = atoms[k]; if (vsite_type[atom] == NOTSET) { - bUsed = FALSE; - for (m = 0; (m < vsnral) && !bUsed; m++) + bool bUsed = false; + for (int m = 0; (m < vsnral) && !bUsed; m++) { GMX_ASSERT(first_atoms != nullptr, "If we've seen a vsite before, we know what its first atom index was"); if (atom == first_atoms[m]) { - bUsed = TRUE; + bUsed = true; bFirstTwo = bFirstTwo && m < 2; } } if (!bUsed) { - bKeep = TRUE; + bKeep = true; } } } @@ -1091,28 +1086,28 @@ static void clean_vsite_bonds(gmx::ArrayRef plist, t_ * a fixed distance due to being constructed from the same * atoms, since this can be numerically unstable. */ - for (m = 0; m < vsnral && !bKeep; m++) /* all constr. atoms */ + for (int m = 0; m < vsnral && !bKeep; m++) /* all constr. atoms */ { at1 = first_atoms[m]; at2 = first_atoms[(m+1) % vsnral]; - bPresent = FALSE; + bool bPresent = false; for (ftype = 0; ftype < F_NRE; ftype++) { if (interaction_function[ftype].flags & IF_CONSTRAINT) { - for (j = 0; (j < plist[ftype].nr) && !bPresent; j++) + for (auto entry = plist[ftype].interactionTypes.begin(); (entry != plist[ftype].interactionTypes.end()) && !bPresent; entry++) { /* all constraints until one matches */ - bPresent = ( ( (plist[ftype].param[j].ai() == at1) && - (plist[ftype].param[j].aj() == at2) ) || - ( (plist[ftype].param[j].ai() == at2) && - (plist[ftype].param[j].aj() == at1) ) ); + bPresent = ( ( (entry->ai() == at1) && + (entry->aj() == at2) ) || + ( (entry->ai() == at2) && + (entry->aj() == at1) ) ); } } } if (!bPresent) { - bKeep = TRUE; + bKeep = true; } } } @@ -1120,32 +1115,30 @@ static void clean_vsite_bonds(gmx::ArrayRef plist, t_ if (bKeep) { - /* now copy the bond to the new array */ - ps->param[kept_i] = ps->param[i]; - kept_i++; + ++bond; } else if (IS_CHEMBOND(cftype)) { - srenew(plist[F_CONNBONDS].param, plist[F_CONNBONDS].nr+1); - plist[F_CONNBONDS].param[plist[F_CONNBONDS].nr] = ps->param[i]; - plist[F_CONNBONDS].nr++; + plist[F_CONNBONDS].interactionTypes.emplace_back(*bond); + bond = ps->interactionTypes.erase(bond); nconverted++; } else { + bond = ps->interactionTypes.erase(bond); nremoved++; } } if (nremoved) { - fprintf(stderr, "Removed %4d %15ss with virtual sites, %5d left\n", - nremoved, interaction_function[cftype].longname, kept_i); + fprintf(stderr, "Removed %4d %15ss with virtual sites, %zu left\n", + nremoved, interaction_function[cftype].longname, ps->size()); } if (nconverted) { - fprintf(stderr, "Converted %4d %15ss with virtual sites to connections, %5d left\n", - nconverted, interaction_function[cftype].longname, kept_i); + fprintf(stderr, "Converted %4d %15ss with virtual sites to connections, %zu left\n", + nconverted, interaction_function[cftype].longname, ps->size()); } if (nOut) { @@ -1157,32 +1150,30 @@ static void clean_vsite_bonds(gmx::ArrayRef plist, t_ nOut, interaction_function[cftype].longname, interaction_function[F_VSITE3OUT].longname ); } - ps->nr = kept_i; } static void clean_vsite_angles(gmx::ArrayRef plist, t_pindex pindex[], int cftype, const int vsite_type[], at2vsitecon_t *at2vc) { - int i, j, k, m, n, nvsite, kept_i; int atom, at1, at2; - bool bKeep, bUsed, bPresent, bAll3FAD, bFirstTwo; InteractionTypeParameters *ps; ps = &(plist[cftype]); - kept_i = 0; - for (i = 0; (i < ps->nr); i++) /* for all angles in the plist */ + int oldSize = ps->size(); + for (auto angle = ps->interactionTypes.begin(); angle != ps->interactionTypes.end(); ) { - int vsnral = 0; - const int *first_atoms = nullptr; + int vsnral = 0; + const int *first_atoms = nullptr; - bKeep = FALSE; - bAll3FAD = TRUE; + bool bKeep = false; + bool bAll3FAD = true; /* check if all virtual sites are constructed from the same atoms */ - nvsite = 0; - for (k = 0; (k < 3) && !bKeep; k++) /* for all atoms in the angle */ + int nvsite = 0; + gmx::ArrayRef atoms = angle->atoms(); + for (int k = 0; (k < 3) && !bKeep; k++) /* for all atoms in the angle */ { - atom = ps->param[i].a[k]; + int atom = atoms[k]; if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN) { nvsite++; @@ -1191,7 +1182,7 @@ static void clean_vsite_angles(gmx::ArrayRef plist, t { /* store construction atoms of first vsite */ vsnral = NRAL(pindex[atom].ftype) - 1; - first_atoms = plist[pindex[atom].ftype].param[pindex[atom].parnr].a + 1; + first_atoms = plist[pindex[atom].ftype].interactionTypes[pindex[atom].parnr].atoms().data() + 1; } else { @@ -1200,28 +1191,28 @@ static void clean_vsite_angles(gmx::ArrayRef plist, t /* check if this vsite is constructed from the same atoms */ if (vsnral == NRAL(pindex[atom].ftype)-1) { - for (m = 0; (m < vsnral) && !bKeep; m++) + for (int m = 0; (m < vsnral) && !bKeep; m++) { - const int *atoms; + const int *subAtoms; - bPresent = FALSE; - atoms = plist[pindex[atom].ftype].param[pindex[atom].parnr].a + 1; - for (n = 0; (n < vsnral) && !bPresent; n++) + bool bPresent = false; + subAtoms = plist[pindex[atom].ftype].interactionTypes[pindex[atom].parnr].atoms().data() + 1; + for (int n = 0; (n < vsnral) && !bPresent; n++) { - if (atoms[m] == first_atoms[n]) + if (subAtoms[m] == first_atoms[n]) { - bPresent = TRUE; + bPresent = true; } } if (!bPresent) { - bKeep = TRUE; + bKeep = true; } } } else { - bKeep = TRUE; + bKeep = true; } } } @@ -1231,30 +1222,30 @@ static void clean_vsite_angles(gmx::ArrayRef plist, t with virtual sites with more than 3 constr. atoms */ if (nvsite == 0 && vsnral > 3) { - bKeep = TRUE; + bKeep = true; } /* check if all non-vsite atoms are used in construction: */ - bFirstTwo = TRUE; - for (k = 0; (k < 3) && !bKeep; k++) /* for all atoms in the angle */ + bool bFirstTwo = true; + for (int k = 0; (k < 3) && !bKeep; k++) /* for all atoms in the angle */ { - atom = ps->param[i].a[k]; + atom = atoms[k]; if (vsite_type[atom] == NOTSET) { - bUsed = FALSE; - for (m = 0; (m < vsnral) && !bUsed; m++) + bool bUsed = false; + for (int m = 0; (m < vsnral) && !bUsed; m++) { GMX_ASSERT(first_atoms != nullptr, "If we've seen a vsite before, we know what its first atom index was"); if (atom == first_atoms[m]) { - bUsed = TRUE; + bUsed = true; bFirstTwo = bFirstTwo && m < 2; } } if (!bUsed) { - bKeep = TRUE; + bKeep = true; } } } @@ -1262,72 +1253,70 @@ static void clean_vsite_angles(gmx::ArrayRef plist, t if (!( bAll3FAD && bFirstTwo ) ) { /* check if all constructing atoms are constrained together */ - for (m = 0; m < vsnral && !bKeep; m++) /* all constr. atoms */ + for (int m = 0; m < vsnral && !bKeep; m++) /* all constr. atoms */ { at1 = first_atoms[m]; at2 = first_atoms[(m+1) % vsnral]; - bPresent = FALSE; - for (j = 0; j < at2vc[at1].nr; j++) + bool bPresent = false; + for (int j = 0; j < at2vc[at1].nr; j++) { if (at2vc[at1].aj[j] == at2) { - bPresent = TRUE; + bPresent = true; } } if (!bPresent) { - bKeep = TRUE; + bKeep = true; } } } if (bKeep) { - /* now copy the angle to the new array */ - ps->param[kept_i] = ps->param[i]; - kept_i++; + ++angle; + } + else + { + angle = ps->interactionTypes.erase(angle); } } - if (ps->nr != kept_i) + if (oldSize != gmx::ssize(*ps)) { - fprintf(stderr, "Removed %4d %15ss with virtual sites, %5d left\n", - ps->nr-kept_i, interaction_function[cftype].longname, kept_i); + fprintf(stderr, "Removed %4zu %15ss with virtual sites, %zu left\n", + oldSize-ps->size(), interaction_function[cftype].longname, ps->size()); } - ps->nr = kept_i; } static void clean_vsite_dihs(gmx::ArrayRef plist, t_pindex pindex[], int cftype, const int vsite_type[]) { - int i, kept_i; InteractionTypeParameters *ps; ps = &(plist[cftype]); - kept_i = 0; - for (i = 0; (i < ps->nr); i++) /* for all dihedrals in the plist */ + int oldSize = ps->size(); + for (auto dih = ps->interactionTypes.begin(); dih != ps->interactionTypes.end(); ) { - int k, m, n, nvsite; - int vsnral = 0; - const int *first_atoms = nullptr; - int atom; - bool bKeep, bUsed, bPresent; - + int vsnral = 0; + const int *first_atoms = nullptr; + int atom; - bKeep = FALSE; + gmx::ArrayRef atoms = dih->atoms(); + bool bKeep = false; /* check if all virtual sites are constructed from the same atoms */ - nvsite = 0; - for (k = 0; (k < 4) && !bKeep; k++) /* for all atoms in the dihedral */ + int nvsite = 0; + for (int k = 0; (k < 4) && !bKeep; k++) /* for all atoms in the dihedral */ { - atom = ps->param[i].a[k]; + atom = atoms[k]; if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN) { if (nvsite == 0) { /* store construction atoms of first vsite */ vsnral = NRAL(pindex[atom].ftype) - 1; - first_atoms = plist[pindex[atom].ftype].param[pindex[atom].parnr].a + 1; + first_atoms = plist[pindex[atom].ftype].interactionTypes[pindex[atom].parnr].atoms().data() + 1; } else { @@ -1336,22 +1325,22 @@ static void clean_vsite_dihs(gmx::ArrayRef plist, t_p /* check if this vsite is constructed from the same atoms */ if (vsnral == NRAL(pindex[atom].ftype)-1) { - for (m = 0; (m < vsnral) && !bKeep; m++) + for (int m = 0; (m < vsnral) && !bKeep; m++) { - const int *atoms; + const int *subAtoms; - bPresent = FALSE; - atoms = plist[pindex[atom].ftype].param[pindex[atom].parnr].a + 1; - for (n = 0; (n < vsnral) && !bPresent; n++) + bool bPresent = false; + subAtoms = plist[pindex[atom].ftype].interactionTypes[pindex[atom].parnr].atoms().data() + 1; + for (int n = 0; (n < vsnral) && !bPresent; n++) { - if (atoms[m] == first_atoms[n]) + if (subAtoms[m] == first_atoms[n]) { - bPresent = TRUE; + bPresent = true; } } if (!bPresent) { - bKeep = TRUE; + bKeep = true; } } } @@ -1366,52 +1355,55 @@ static void clean_vsite_dihs(gmx::ArrayRef plist, t_p /* keep all dihedrals with no virtual sites in them */ if (nvsite == 0) { - bKeep = TRUE; + bKeep = true; } /* check if all atoms in dihedral are either virtual sites, or used in construction of virtual sites. If so, keep it, if not throw away: */ - for (k = 0; (k < 4) && !bKeep; k++) /* for all atoms in the dihedral */ + for (int k = 0; (k < 4) && !bKeep; k++) /* for all atoms in the dihedral */ { GMX_ASSERT(vsnral != 0, "If we've seen a vsite before, we know how many constructing atoms it had"); GMX_ASSERT(first_atoms != nullptr, "If we've seen a vsite before, we know what its first atom index was"); - atom = ps->param[i].a[k]; + atom = atoms[k]; if (vsite_type[atom] == NOTSET) { /* vsnral will be set here, we don't get here with nvsite==0 */ - bUsed = FALSE; - for (m = 0; (m < vsnral) && !bUsed; m++) + bool bUsed = false; + for (int m = 0; (m < vsnral) && !bUsed; m++) { if (atom == first_atoms[m]) { - bUsed = TRUE; + bUsed = true; } } if (!bUsed) { - bKeep = TRUE; + bKeep = true; } } } if (bKeep) { - ps->param[kept_i] = ps->param[i]; - kept_i++; + ++dih; } + else + { + dih = ps->interactionTypes.erase(dih); + } + } - if (ps->nr != kept_i) + if (oldSize != gmx::ssize(*ps)) { - fprintf(stderr, "Removed %4d %15ss with virtual sites, %5d left\n", - ps->nr-kept_i, interaction_function[cftype].longname, kept_i); + fprintf(stderr, "Removed %4zu %15ss with virtual sites, %zu left\n", + oldSize-ps->size(), interaction_function[cftype].longname, ps->size()); } - ps->nr = kept_i; } void clean_vsite_bondeds(gmx::ArrayRef plist, int natoms, bool bRmVSiteBds) { - int i, k, nvsite, ftype, vsite, parnr; + int nvsite, vsite; int *vsite_type; t_pindex *pindex; at2vsitecon_t *at2vc; @@ -1419,20 +1411,20 @@ void clean_vsite_bondeds(gmx::ArrayRef plist, int nat pindex = nullptr; /* avoid warnings */ /* make vsite_type array */ snew(vsite_type, natoms); - for (i = 0; i < natoms; i++) + for (int i = 0; i < natoms; i++) { vsite_type[i] = NOTSET; } nvsite = 0; - for (ftype = 0; ftype < F_NRE; ftype++) + for (int ftype = 0; ftype < F_NRE; ftype++) { if (interaction_function[ftype].flags & IF_VSITE) { - nvsite += plist[ftype].nr; - i = 0; - while (i < plist[ftype].nr) + nvsite += plist[ftype].size(); + int i = 0; + while (i < gmx::ssize(plist[ftype])) { - vsite = plist[ftype].param[i].ai(); + vsite = plist[ftype].interactionTypes[i].ai(); if (vsite_type[vsite] == NOTSET) { vsite_type[vsite] = ftype; @@ -1443,7 +1435,7 @@ void clean_vsite_bondeds(gmx::ArrayRef plist, int nat } if (ftype == F_VSITEN) { - while (i < plist[ftype].nr && plist[ftype].param[i].ai() == vsite) + while (i < gmx::ssize(plist[ftype]) && plist[ftype].interactionTypes[i].ai() == vsite) { i++; } @@ -1466,7 +1458,7 @@ void clean_vsite_bondeds(gmx::ArrayRef plist, int nat at2vc = make_at2vsitecon(natoms, plist); snew(pindex, natoms); - for (ftype = 0; ftype < F_NRE; ftype++) + for (int ftype = 0; ftype < F_NRE; ftype++) { /* Here we skip VSITEN. In neary all practical use cases this * is not an issue, since VSITEN is intended for constructing @@ -1482,9 +1474,9 @@ void clean_vsite_bondeds(gmx::ArrayRef plist, int nat if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN) { - for (parnr = 0; (parnr < plist[ftype].nr); parnr++) + for (int parnr = 0; (parnr < gmx::ssize(plist[ftype])); parnr++) { - k = plist[ftype].param[parnr].ai(); + int k = plist[ftype].interactionTypes[parnr].ai(); pindex[k].ftype = ftype; pindex[k].parnr = parnr; } @@ -1492,7 +1484,7 @@ void clean_vsite_bondeds(gmx::ArrayRef plist, int nat } /* remove interactions that include virtual sites */ - for (ftype = 0; ftype < F_NRE; ftype++) + for (int ftype = 0; ftype < F_NRE; ftype++) { if ( ( ( interaction_function[ftype].flags & IF_BOND ) && bRmVSiteBds ) || ( interaction_function[ftype].flags & IF_CONSTRAINT ) ) @@ -1512,7 +1504,7 @@ void clean_vsite_bondeds(gmx::ArrayRef plist, int nat } } /* check that no remaining constraints include virtual sites */ - for (ftype = 0; ftype < F_NRE; ftype++) + for (int ftype = 0; ftype < F_NRE; ftype++) { if (interaction_function[ftype].flags & IF_CONSTRAINT) { diff --git a/src/gromacs/gmxpreprocess/x2top.cpp b/src/gromacs/gmxpreprocess/x2top.cpp index 3a587ab5c3..6c54e707c8 100644 --- a/src/gromacs/gmxpreprocess/x2top.cpp +++ b/src/gromacs/gmxpreprocess/x2top.cpp @@ -97,21 +97,12 @@ static void mk_bonds(int nnm, t_nm2type nmt[], t_atoms *atoms, const rvec x[], InteractionTypeParameters *bond, int nbond[], bool bPBC, matrix box) { - t_param b; - int i, j; - t_pbc pbc; - rvec dx; - real dx2; - - for (i = 0; (i < MAXATOMLIST); i++) - { - b.a[i] = -1; - } - for (i = 0; (i < MAXFORCEPARAM); i++) - { - b.c[i] = 0.0; - } + int i, j; + t_pbc pbc; + rvec dx; + real dx2; + std::array forceParam = {0.0}; if (bPBC) { set_pbc(&pbc, -1, box); @@ -138,10 +129,9 @@ static void mk_bonds(int nnm, t_nm2type nmt[], if (is_bond(nnm, nmt, *atoms->atomname[i], *atoms->atomname[j], std::sqrt(dx2))) { - b.ai() = i; - b.aj() = j; - b.c0() = std::sqrt(dx2); - add_param_to_list (bond, &b); + forceParam[0] = std::sqrt(dx2); + std::vector atoms = {i, j}; + add_param_to_list (bond, InteractionType(atoms, forceParam)); nbond[i]++; nbond[j]++; } @@ -206,7 +196,7 @@ static void lo_set_force_const(InteractionTypeParameters *plist, real c[], int n double cc; char buf[32]; - for (int i = 0; (i < plist->nr); i++) + for (auto ¶m : plist->interactionTypes) { if (!bParam) { @@ -219,13 +209,13 @@ static void lo_set_force_const(InteractionTypeParameters *plist, real c[], int n { if (bRound) { - sprintf(buf, "%.2e", plist->param[i].c[0]); + sprintf(buf, "%.2e", param.c0()); sscanf(buf, "%lf", &cc); c[0] = cc; } else { - c[0] = plist->param[i].c[0]; + c[0] = param.c0(); } if (bDih) { @@ -240,12 +230,13 @@ static void lo_set_force_const(InteractionTypeParameters *plist, real c[], int n } } GMX_ASSERT(nrfp <= MAXFORCEPARAM/2, "Only 6 parameters may be used for an interaction"); + std::array forceParam; for (int j = 0; (j < nrfp); j++) { - plist->param[i].c[j] = c[j]; - plist->param[i].c[nrfp+j] = c[j]; + forceParam[j] = c[j]; + forceParam[nrfp+j] = c[j]; } - set_p_string(&(plist->param[i]), ""); + param = InteractionType(param.atoms(), forceParam); } } @@ -276,24 +267,24 @@ static void calc_angles_dihs(InteractionTypeParameters *ang, InteractionTypePara { set_pbc(&pbc, epbcXYZ, box); } - for (int i = 0; (i < ang->nr); i++) + for (auto &angle : ang->interactionTypes) { - int ai = ang->param[i].ai(); - int aj = ang->param[i].aj(); - int ak = ang->param[i].ak(); + int ai = angle.ai(); + int aj = angle.aj(); + int ak = angle.ak(); real th = RAD2DEG*bond_angle(x[ai], x[aj], x[ak], bPBC ? &pbc : nullptr, r_ij, r_kj, &costh, &t1, &t2); - ang->param[i].c0() = th; + angle.setForceParameter(0, th); } - for (int i = 0; (i < dih->nr); i++) + for (auto dihedral : dih->interactionTypes) { - int ai = dih->param[i].ai(); - int aj = dih->param[i].aj(); - int ak = dih->param[i].ak(); - int al = dih->param[i].al(); + int ai = dihedral.ai(); + int aj = dihedral.aj(); + int ak = dihedral.ak(); + int al = dihedral.al(); real ph = RAD2DEG*dih_angle(x[ai], x[aj], x[ak], x[al], bPBC ? &pbc : nullptr, r_ij, r_kj, r_kl, m, n, &t1, &t2, &t3); - dih->param[i].c0() = ph; + dihedral.setForceParameter(0, ph); } } @@ -308,23 +299,24 @@ static void dump_hybridization(FILE *fp, t_atoms *atoms, int nbonds[]) static void print_pl(FILE *fp, gmx::ArrayRef plist, int ftp, const char *name, char ***atomname) { - if (plist[ftp].nr > 0) + if (!plist[ftp].interactionTypes.empty()) { fprintf(fp, "\n"); fprintf(fp, "[ %s ]\n", name); - int nral = interaction_function[ftp].nratoms; int nrfp = interaction_function[ftp].nrfpA; - for (int i = 0; (i < plist[ftp].nr); i++) + for (const auto ¶m : plist[ftp].interactionTypes) { - for (int j = 0; (j < nral); j++) + gmx::ArrayRef atoms = param.atoms(); + gmx::ArrayRef forceParam = param.forceParam(); + for (const auto &atom : atoms) { - fprintf(fp, " %5s", *atomname[plist[ftp].param[i].a[j]]); + fprintf(fp, " %5s", *atomname[atom]); } for (int j = 0; (j < nrfp); j++) { - if (plist[ftp].param[i].c[j] != NOTSET) + if (forceParam[j] != NOTSET) { - fprintf(fp, " %10.3e", plist[ftp].param[i].c[j]); + fprintf(fp, " %10.3e", forceParam[j]); } } fprintf(fp, "\n"); @@ -542,15 +534,15 @@ int gmx_x2top(int argc, char *argv[]) if (!bPairs) { - plist[F_LJ14].nr = 0; + plist[F_LJ14].interactionTypes.clear(); } fprintf(stderr, - "There are %4d %s dihedrals, %4d impropers, %4d angles\n" - " %4d pairs, %4d bonds and %4d atoms\n", - plist[F_PDIHS].nr, + "There are %4zu %s dihedrals, %4zu impropers, %4zu angles\n" + " %4zu pairs, %4zu bonds and %4d atoms\n", + plist[F_PDIHS].size(), bOPLS ? "Ryckaert-Bellemans" : "proper", - plist[F_IDIHS].nr, plist[F_ANGLES].nr, - plist[F_LJ14].nr, plist[F_BONDS].nr, atoms->nr); + plist[F_IDIHS].size(), plist[F_ANGLES].size(), + plist[F_LJ14].size(), plist[F_BONDS].size(), atoms->nr); calc_angles_dihs(&plist[F_ANGLES], &plist[F_PDIHS], x, bPBC, box); -- 2.22.0