#include "hackblock.h"
-static void clear_atom_list(int i0, int a[])
+void add_param(InteractionTypeParameters *ps,
+ int ai,
+ int aj,
+ gmx::ArrayRef<const real> 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<int> atoms = {ai, aj};
+ std::vector<real> 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<int> atoms = {ai, aj, ak, al};
+ std::vector<real> 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<int> atoms = {ai, aj, ak, al};
+ std::vector<real> 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<int> 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<int> 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<int> atoms = {ai, aj, ak};
+ std::vector<real> 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<int> atoms = {ai, aj, ak, al};
+ std::vector<real> 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<int> 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<int> atoms = {ai, aj, ak, al, am};
+ ps->interactionTypes.emplace_back(InteractionType(atoms, {}));
}
int search_jtype(const PreprocessResidue &localPpResidue, const char *name, bool bNterm)
#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<const real> 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);
*/
static int
assign_param(t_functype ftype, t_iparams *newparam,
- real old[MAXFORCEPARAM], int comb, double reppow)
+ gmx::ArrayRef<const real> old, int comb, double reppow)
{
bool all_param_zero = true;
}
static int enter_params(gmx_ffparams_t *ffparams, t_functype ftype,
- real forceparams[MAXFORCEPARAM], int comb, real reppow,
+ gmx::ArrayRef<const real> forceparams, int comb, real reppow,
int start, bool bAppend)
{
t_iparams newparam;
}
static void append_interaction(InteractionList *ilist,
- int type, int nral, const int a[MAXATOMLIST])
+ int type, gmx::ArrayRef<const int> 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);
}
}
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());
}
}
}
gmx::ArrayRef<const InteractionTypeParameters> plist = intermolecular_interactions->plist;
- if (plist[i].nr > 0)
+ if (!plist[i].interactionTypes.empty())
{
flags = interaction_function[i].flags;
/* For intermolecular interactions we (currently)
#include <cstring>
#include <algorithm>
+#include <numeric>
#include "gromacs/fileio/confio.h"
#include "gromacs/gmxpreprocess/gpp_nextnb.h"
#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<const real> 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<const t_param *>(a1);
- p2 = static_cast<const t_param *>(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<const t_param *>(a1);
- p2 = static_cast<const t_param *>(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<const t_param *>(d1);
- p2 = static_cast<const t_param *>(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<InteractionType> *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<InteractionTypeParameters> plist, int ftype)
+static void cppar(gmx::ArrayRef<const InteractionType> types,
+ gmx::ArrayRef<InteractionTypeParameters> 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<const t_param *>(a);
- pb = static_cast<const t_param *>(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<InteractionType> 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<size_t>(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<const int> 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<int>(strlen(aname)) > 1) && (c0 >= '0') && (c0 <= '9'))
{
- c1 = toupper(aname[1]);
+ char c1 = toupper(aname[1]);
if (c1 == 'H')
{
nh++;
/* 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<InteractionType> clean_dih(gmx::ArrayRef<const InteractionType> dih,
+ gmx::ArrayRef<const InteractionType> 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<InteractionType, int>(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<InteractionType, int>(*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);
}
}
* 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;
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<InteractionType> 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<MoleculePatchDatabase> globalPatches, t_param **improper,
- bool bAllowMissing)
+static std::vector<InteractionType> get_impropers(t_atoms *atoms,
+ gmx::ArrayRef<MoleculePatchDatabase> globalPatches,
+ bool bAllowMissing)
{
- int nimproper, start, ninc, nalloc;
- int ai[MAXATOMLIST];
-
- ninc = 500;
- nalloc = ninc;
- snew(*improper, nalloc);
+ std::vector<InteractionType> 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++)
BondedInteractionList *impropers = &globalPatches[i].rb[ebtsIDIHS];
for (const auto &bondeds : impropers->b)
{
- bool bStop = false;
+ bool bStop = false;
+ std::vector<int> 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))
}
}
- return nimproper;
+ return improper;
}
static int nb_dist(t_nextnb *nnb, int ai, int aj)
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<std::string> anm,
+ int resind, t_atoms *atoms, gmx::ArrayRef<const int> 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]]));
}
}
gmx::ArrayRef<InteractionTypeParameters> plist, t_excls excls[], gmx::ArrayRef<MoleculePatchDatabase> 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<InteractionType> ang;
+ std::vector<InteractionType> dih;
+ std::vector<InteractionType> pai;
+ std::vector<InteractionType> improper;
+
+ std::array<std::string, 4> anm;
if (!globalPatches.empty())
{
/* Generate every angle only once */
if (i < k1)
{
- if (nang == maxang)
- {
- maxang += ninc;
- srenew(ang, maxang);
- }
- ang[nang].ai() = i;
- ang[nang].aj() = j1;
- ang[nang].ak() = k1;
- ang[nang].c0() = NOTSET;
- ang[nang].c1() = NOTSET;
- set_p_string(&(ang[nang]), "");
+ std::vector<int> 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)
{
}
if (bFound)
{
- set_p_string(&(ang[nang]), bondeds.s.c_str());
+ name = bondeds.s;
/* Mark that we found a match for this entry */
bondeds.match = true;
}
}
while (res < maxres);
}
- nang++;
+ ang.push_back(InteractionType(atomNumbers, {}, name));
}
/* Generate every dihedral, 1-4 exclusion and 1-4 interaction
only once */
int l1 = nnb->a[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<int> 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)
{
}
if (bFound)
{
- set_p_string(&dih[ndih], bondeds.s.c_str());
+ name = bondeds.s;
/* Mark that we found a match for this entry */
bondeds.match = true;
/* Set the last parameter to be able to see
if the dihedral was in the rtp list.
*/
- dih[ndih].c[MAXFORCEPARAM-1] = DIHEDRAL_WAS_SET_IN_RTP;
nFound++;
- ndih++;
- /* Set the next direct in case the rtp contains
- multiple entries for this dihedral.
- */
- 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;
- }
+ dih.push_back(InteractionType(atomNumbers, {}, name));
+ dih.back().setForceParameter(MAXFORCEPARAM-1, DIHEDRAL_WAS_SET_IN_RTP);
}
}
}
}
if (nFound == 0)
{
- 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]), "");
- ndih++;
+ std::vector<int> atoms = {i, j1, k1, l1};
+ dih.push_back(InteractionType(atoms, {}, ""));
}
int nbd = nb_dist(nnb, i, l1);
if (rtpFFDB[0].bGenerateHH14Interactions ||
!(is_hydro(atoms, i1) && is_hydro(atoms, i2)))
{
- if (npai == maxpai)
- {
- maxpai += ninc;
- srenew(pai, maxpai);
- }
- pai[npai].ai() = i1;
- pai[npai].aj() = i2;
- pai[npai].c0() = NOTSET;
- pai[npai].c1() = NOTSET;
- set_p_string(&(pai[npai]), "");
- npai++;
+ std::vector<int> atoms = {i1, i2};
+ pai.push_back(InteractionType(atoms, {}, ""));
}
}
}
continue;
}
/* Hm - entry not used, let's see if we can find all atoms */
- if (nang == maxang)
- {
- maxang += ninc;
- srenew(ang, maxang);
- }
- bool bFound = true;
+ std::vector<int> atomNumbers;
+ bool bFound = true;
for (int k = 0; k < 3 && bFound; k++)
{
const char *p = bondeds.a[k].c_str();
p++;
res++;
}
- ang[nang].a[k] = search_res_atom(p, res, atoms, "angle", TRUE);
- bFound = (ang[nang].a[k] != -1);
+ atomNumbers.emplace_back(search_res_atom(p, res, atoms, "angle", TRUE));
+ bFound = (atomNumbers.back() != -1);
}
- ang[nang].c0() = NOTSET;
- ang[nang].c1() = NOTSET;
if (bFound)
{
- set_p_string(&(ang[nang]), bondeds.s.c_str());
bondeds.match = true;
/* Incrementing nang means we save this angle */
- nang++;
+ ang.push_back(InteractionType(atomNumbers, {}, bondeds.s));
}
}
continue;
}
/* Hm - entry not used, let's see if we can find all atoms */
- if (ndih == maxdih)
- {
- maxdih += ninc;
- srenew(dih, maxdih);
- }
- bool bFound = true;
+ std::vector<int> atomNumbers;
+ bool bFound = true;
for (int k = 0; k < 4 && bFound; k++)
{
const char *p = bondeds.a[k].c_str();
p++;
res++;
}
- dih[ndih].a[k] = search_res_atom(p, res, atoms, "dihedral", TRUE);
- bFound = (dih[ndih].a[k] != -1);
- }
- for (int m = 0; m < MAXFORCEPARAM; m++)
- {
- dih[ndih].c[m] = NOTSET;
+ atomNumbers.emplace_back(search_res_atom(p, res, atoms, "dihedral", TRUE));
+ bFound = (atomNumbers.back() != -1);
}
if (bFound)
{
- set_p_string(&(dih[ndih]), bondeds.s.c_str());
bondeds.match = true;
/* Incrementing ndih means we save this dihedral */
- ndih++;
+ dih.push_back(InteractionType(atomNumbers, {}, bondeds.s));
}
}
}
}
/* Sort angles with respect to j-i-k (middle atom first) */
- if (nang > 1)
+ if (ang.size() > 1)
{
- qsort(ang, nang, static_cast<size_t>(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<size_t>(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<size_t>(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);
}
/* 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)
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)
{
{
/* 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) )
{
char tpname[32], nexttpname[32];
int *o2n, *newvsite_type, *newcgnr, ats[MAXATOMSPERRESIDUE];
t_atom *newatom;
- InteractionTypeParameters *params;
char ***newatomname;
char *resnm = nullptr;
int cmplength;
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<const int> atomNumbers(parm->atoms());
+ std::vector<int> 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 */
/* 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,
{
/* 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)
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
{
return NOTSET;
}
+ gmx::ArrayRef<const real> 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()
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)
}
}
-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))
{
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;
fprintf (out, "\n");
}
-static int search_atomtypes(const PreprocessingAtomTypes *ga,
- int *n,
- gmx::ArrayRef<int> typelist,
- int thistype,
- t_param param[],
- int ftype)
+static int search_atomtypes(const PreprocessingAtomTypes *ga,
+ int *n,
+ gmx::ArrayRef<int> typelist,
+ int thistype,
+ gmx::ArrayRef<const InteractionType> 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)
for (int j = 0; j < ntype && bFound; j++)
{
/* Check nonbonded parameters */
- for (int k = 0; k < nrfp && bFound; k++)
+ gmx::ArrayRef<const real> forceParam1 = interactionTypes[ntype*typelist[i]+j].forceParam();
+ gmx::ArrayRef<const real> 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 */
*/
/* Get nonbonded interaction type */
- if (plist[F_LJ].nr > 0)
+ if (plist[F_LJ].size() > 0)
{
ftype = F_LJ;
}
{
atoms->atom[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);
}
}
if (wall_atomtype[i] >= 0)
{
wall_atomtype[i] = search_atomtypes(this, &nat, typelist, wall_atomtype[i],
- plist[ftype].param, ftype);
+ plist[ftype].interactionTypes, ftype);
}
}
/* We now have a list of unique atomtypes in typelist */
/* Renumber nlist */
- /* Renumber nlist */
- t_param *nbsnew = nullptr;
- snew(nbsnew, plist[ftype].nr);
+ std::vector<InteractionType> 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
struct gmx_mtop_t;
struct t_atom;
struct t_atomtypes;
-struct t_param;
+class InteractionType;
struct InteractionTypeParameters;
struct t_symtab;
* \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.
* \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.
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",
s[(*nrf)++].aj = aj;
s[(*nrf)].aj = ai;
s[(*nrf)++].ai = aj;
+ i++;
}
}
if (IS_CHEMBOND(i))
{
/* we need every bond twice (bidirectional) */
- nrbonds += 2*plist[i].nr;
+ nrbonds += 2*plist[i].size();
}
}
#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<const int> atoms,
+ gmx::ArrayRef<const real> 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);
void MoleculeInformation::partialCleanUp()
{
done_block(&mols);
- done_plist(plist);
}
void MoleculeInformation::fullCleanUp()
done_atom (&atoms);
done_block(&cgs);
done_block(&mols);
- done_plist(plist);
}
static int rm_interactions(int ifunc, gmx::ArrayRef<MoleculeInformation> mols)
/* 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;
}
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;
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);
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;
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",
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];
}
}
}
/* 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",
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];
}
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]);
}
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]);
}
{
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;
}
if (!molb.posres_xA.empty() || !molb.posres_xB.empty())
{
std::vector<gmx::RVec> &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];
}
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];
}
gmx::ArrayRef<const MoleculeInformation> mi,
warninp *wi)
{
- int count, count_mol, i;
+ int count, count_mol;
char buf[STRLEN];
count = 0;
count_mol = 0;
gmx::ArrayRef<const InteractionTypeParameters> 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();
}
}
}
std::array<InteractionTypeParameters, F_NRE> plist;
- init_plist(plist);
gmx_mtop_t sys;
PreprocessingAtomTypes atypes;
if (debug)
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
#include <string>
+#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<const int> atoms,
+ gmx::ArrayRef<const real> 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<int> atoms() { return atoms_; }
+ //! Const view on all atoms numbers that are actually set.
+ gmx::ArrayRef<const int> atoms() const { return atoms_; }
+ //! View on all of the force field parameters
+ gmx::ArrayRef<const real> forceParam() const { return forceParam_; }
+ //! View on all of the force field parameters
+ gmx::ArrayRef<real> 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<int> atoms_;
+ //! Force parameters (eg. b0 = forceParam[0])
+ std::array<real, MAXFORCEPARAM> 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<InteractionType> 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<real> cmap;
+ std::vector<real> cmap;
//! The five atomtypes followed by a number that identifies the type.
- std::vector<int> cmapAtomTypes;
+ std::vector<int> 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
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++)
{
{
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;
{
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;
sfree(bbb);
sfree(n_mask);
sfree(m_mask);
- sfree(param);
return nresolved;
}
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);
}
}
fprintf(stderr, "Warning: Short Bond (%d-%d = %g nm)\n",
ai+1, aj+1, std::sqrt(dist2));
}
- add_param(psb, ai, aj, nullptr, patch.s.c_str());
+ add_param(psb, ai, aj, {}, patch.s.c_str());
}
}
/* add bonds from list of hacks (each added atom gets a bond) */
{
switch (patch.tp)
{
- case 9 : /* COOH terminus */
- add_param(psb, i, i+1, nullptr, nullptr); /* C-O */
- add_param(psb, i, i+2, nullptr, nullptr); /* C-OA */
- add_param(psb, i+2, i+3, nullptr, nullptr); /* OA-H */
+ case 9 : /* COOH terminus */
+ add_param(psb, i, i+1, {}, nullptr); /* C-O */
+ add_param(psb, i, i+2, {}, nullptr); /* C-OA */
+ add_param(psb, i+2, i+3, {}, nullptr); /* OA-H */
break;
default:
for (int k = 0; (k < patch.nr); k++)
{
- add_param(psb, i, i+k+1, nullptr, nullptr);
+ add_param(psb, i, i+k+1, {}, nullptr);
}
}
}
}
}
-static int pcompar(const void *a, const void *b)
+static bool pcompar(const InteractionType &a, const InteractionType &b)
{
- const t_param *pa, *pb;
- int d;
- pa = static_cast<const t_param *>(a);
- pb = static_cast<const t_param *>(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<size_t>(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
{
int i, nmissat;
int bts[ebtsNR];
- init_plist(plist);
ResidueType rt;
/* Make bonds */
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());
}
}
/* 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);
/* 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);
double m;
int nratt = 0;
t_atom *a;
- t_param *nb;
std::vector<std::string> files = fflib_search_file_end(ffdir, ".atp", TRUE);
snew(a, 1);
- snew(nb, 1);
- PreprocessingAtomTypes at;
+ PreprocessingAtomTypes at;
for (const auto &filename : files)
{
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);
}
int n2m;
t_2morse *t2m;
- bool *bRemoveHarm;
-
/* First get the data */
t2m = read_dissociation_energies(&n2m);
if (n2m <= 0)
/* 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<int> atoms = {ni, nj};
+ std::vector<real> 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++;
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<int>(std::sqrt(static_cast<double>(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))
{
}
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<int, 2> atomNumbers;
+ std::array<real, MAXFORCEPARAM> 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<const real> 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.
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++;
}
}
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<real, MAXFORCEPARAM> forceParam = {NOTSET};
/* Fill the matrix with force parameters */
switch (ftype)
{
{
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));
}
}
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);
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));
}
}
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++;
}
const char *entry;
int ptype;
} t_xlate;
- t_xlate xl[eptNR] = {
+ t_xlate xl[eptNR] = {
{ "A", eptAtom },
{ "N", eptNucleus },
{ "S", eptShell },
{ "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",
{
c[j] = 0.0;
}
+ std::array<real, MAXFORCEPARAM> forceParam;
if (strlen(type) == 1 && isdigit(type[0]))
{
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);
{
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);
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.
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
if (bAllowRepeat && nr > 1)
{
isContinuationOfBlock = true;
+ gmx::ArrayRef<const int> newParAtom = b.atoms();
+ gmx::ArrayRef<const int> 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;
}
bool haveErrored = false;
for (int i = 0; (i < nr); i++)
{
- gmx::ArrayRef<const int> bParams(b->a, b->a + nral);
- gmx::ArrayRef<const int> testParams(bt->param[i].a, bt->param[i].a + nral);
+ gmx::ArrayRef<const int> bParams = b.atoms();
+ gmx::ArrayRef<const int> 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)
{
warning(wi, message);
fprintf(stderr, " old: ");
- for (int j = 0; (j < nrfp); j++)
+ gmx::ArrayRef<const real> 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);
/* 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<const real> 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]);
}
}
}
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<real> 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
*/
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<int> atoms;
+ gmx::ArrayRef<const int> 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()));
}
}
char *line,
warninp *wi)
{
- const char *formal[MAXATOMLIST+1] = {
+ const char *formal[MAXATOMLIST+1] = {
"%s",
"%s%s",
"%s%s%s",
"%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",
"%*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))
{
}
}
}
- for (i = 0; (i < nral); i++)
+ std::vector<int> atoms;
+ std::array<real, MAXFORCEPARAM> 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);
}
gpp_bond_atomtype *bat, char *line,
warninp *wi)
{
- const char *formal[MAXATOMLIST+1] = {
+ const char *formal[MAXATOMLIST+1] = {
"%s",
"%s%s",
"%s%s%s",
"%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",
"%*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",
"%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.
*
}
}
- for (i = 0; (i < 4); i++)
+ std::vector<int> atoms;
+ std::array<real, MAXFORCEPARAM> 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);
}
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;
bt[F_CMAP].cmapAngles++; /* Since we are incrementing here, we need to subtract later, see (*****) */
nct = (nral+1) * bt[F_CMAP].cmapAngles;
+ std::vector<int> 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));
warning_error(wi, message);
}
- /* Is this correct?? */
- for (int i = 0; i < MAXFORCEPARAM; i++)
- {
- p.c[i] = NOTSET;
- }
+ std::array<real, MAXFORCEPARAM> 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);
}
mol->back().excl_set = false;
}
+static bool findIfAllNBAtomsMatch(gmx::ArrayRef<const int> atomsFromParameterArray,
+ gmx::ArrayRef<const int> 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<InteractionTypeParameters> 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))
{
GMX_ASSERT(ntype * ntype == nr, "Number of pairs of generated non-bonded parameters should be a perfect square");
if (bB)
{
- ti = at->atom[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<const int> 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<const real> 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;
static bool default_cmap_params(gmx::ArrayRef<InteractionTypeParameters> bondtype,
t_atoms *at, PreprocessingAtomTypes *atypes,
- t_param *p, bool bB,
+ InteractionType *p, bool bB,
int *cmap_type, int *nparam_def,
warninp *wi)
{
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;
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);
}
/* 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
{
}
}
-static bool defaultInteractionTypeParameters(int ftype, gmx::ArrayRef<InteractionTypeParameters> 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<const int> atomsFromParameterArray,
+ gmx::ArrayRef<const int> 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<InteractionType>::iterator
+defaultInteractionTypeParameters(int ftype, gmx::ArrayRef<InteractionTypeParameters> 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
* The rule in that case is that additional matches
* HAVE to be on adjacent lines!
*/
- bSame = TRUE;
- /* Continue from current i value */
- for (j = i + 2; j < nr && bSame; j += 2)
+ bool bSame = true;
+ /* Continue from current iterator position */
+ for (auto nextPos = prevPos + 2; (nextPos < bt[ftype].interactionTypes.end()) && bSame; nextPos += 2)
{
- pj = &(bt[ftype].param[j]);
- bSame = (pi->ai() == 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++;
/* nparam_found will be increased as long as the numbers match */
}
}
+ *nparam_def = nparam_found;
+ return prevPos;
}
else /* Not a dihedral */
{
- int i, j;
-
- for (i = 0; ((i < nr) && !bFound); i++)
- {
- pi = &(bt[ftype].param[i]);
- if (bB)
- {
- for (j = 0; ((j < nral) &&
- (atypes->bondAtomTypeFromAtomType(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<const int> 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;
}
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",
"%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",
"%*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;
}
/* 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)
{
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])
}
/* default force parameters */
- for (j = 0; (j < MAXATOMLIST); j++)
+ std::vector<int> 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<real, MAXFORCEPARAM> forceParam = {0.0};
/* Get force params for normal and free energy perturbation
* studies, as determined by types!
*/
+ InteractionType param(atoms, forceParam, "");
+ std::vector<InteractionType>::iterator foundAParameter = bondtype[ftype].interactionTypes.end();
+ std::vector<InteractionType>::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<const real> 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<const real> 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)
}
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;
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<const int> 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)
/* 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];
}
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))
{
/* 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"
{
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
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;
}
}
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<const int> 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)
}
}
+ gmx::ArrayRef<const real> 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++;
}
}
/* 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,
*/
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<const real> 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);
}
}
}
"%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);
}
/* 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)
{
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])
{
}
/* default force parameters */
- for (j = 0; (j < MAXATOMLIST); j++)
- {
- param.a[j] = aa[j]-1;
- }
- for (j = 0; (j < MAXFORCEPARAM); j++)
+ std::vector<int> atoms;
+ for (int j = 0; (j < nral); j++)
{
- param.c[j] = 0.0;
+ atoms.emplace_back(aa[j]-1);
}
-
+ std::array<real, MAXFORCEPARAM> 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);
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);
}
}
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<real, MAXFORCEPARAM> forceParam = {0.0};
ptr = line;
ret = sscanf(ptr, "%d%n", &a, &n);
ptr += n;
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;
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<int> 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);
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;
* 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<real, MAXFORCEPARAM> 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);
static void convert_pairs_to_pairsQ(gmx::ArrayRef<InteractionTypeParameters> 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<InteractionType> paramnew;
- paramp1 = plist[F_LJ14].param;
- paramp2 = plist[F_LJC14_Q].param;
+ gmx::ArrayRef<const InteractionType> paramp1 = plist[F_LJ14].interactionTypes;
+ gmx::ArrayRef<const InteractionType> 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
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<real> 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<int>(std::sqrt(static_cast<double>(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<int>(std::sqrt(static_cast<double>(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)
{
"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<int> atoms = {i, j};
+ std::vector<real> 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));
}
}
}
struct t_block;
struct MoleculeInformation;
struct t_nbparam;
-struct t_param;
+class InteractionType;
struct InteractionTypeParameters;
struct PreprocessResidue;
struct warninp;
#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<const int> a)
{
int nh;
{
char ***info = atoms->atomname;
real b_ij, b_jk;
- int i, j;
-
if (nshake != eshNONE)
{
switch (nshake)
{
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<int> 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<real> 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 */
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<int> atomNumbers = {parm->ai(), parm->aj()};
+ std::vector<real> 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;
}
}
}
/* 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<int>(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<InteractionTypeParameters> 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<InteractionTypeParameters> 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 */
const InteractionTypeParameters *bt = &(plist[ftype]);
- if (bt->nr == 0)
+ if (bt->size() == 0)
{
return;
}
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<const int> 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<const real> 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]);
}
}
int ftype, int fsubtype, gmx::ArrayRef<const InteractionTypeParameters> 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);
}
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<InteractionTypeParameters> plist);
-void done_plist(gmx::ArrayRef<InteractionTypeParameters> plist);
/* PRINTING */
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<vsitebondparam_t> vSiteBondedParameters;
+ std::vector<VsiteBondParameter> vSiteBondedParameters;
};
typedef struct {
}
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<const int> 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<const int> atoms,
gmx::ArrayRef<const Atom2VsiteBond> at2vb,
int *nrbond, t_mybonded **bonds,
int *nrang, t_mybonded **angles,
{
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;
}
}
}
make_at2vsitebond(int natoms, gmx::ArrayRef<InteractionTypeParameters> plist)
{
bool *bVSI;
- t_iatom *aa;
std::vector<Atom2VsiteBond> at2vb(natoms);
{
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<const int> 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;
}
}
}
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<const int> 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]);
}
}
}
static at2vsitecon_t *make_at2vsitecon(int natoms, gmx::ArrayRef<InteractionTypeParameters> 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<const int> 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);
static void done_at2vsitecon(int natoms, at2vsitecon_t *at2vc)
{
- int i;
-
- for (i = 0; i < natoms; i++)
+ for (int i = 0; i < natoms; i++)
{
if (at2vc[i].nr)
{
int nrang, t_mybonded *angles,
int nridih, t_mybonded *idihs )
{
- int i;
-
if (nrbond)
{
fprintf(fp, "bonds:");
- for (i = 0; i < nrbond; i++)
+ for (int i = 0; i < nrbond; i++)
{
fprintf(fp, " %d-%d (%g)",
bonds[i].ai()+1, bonds[i].aj()+1, bonds[i].c);
if (nrang)
{
fprintf(fp, "angles:");
- for (i = 0; i < nrang; i++)
+ for (int i = 0; i < nrang; i++)
{
fprintf(fp, " %d-%d-%d (%g)",
angles[i].ai()+1, angles[i].aj()+1,
if (nridih)
{
fprintf(fp, "idihs:");
- for (i = 0; i < nridih; i++)
+ for (int i = 0; i < nridih; i++)
{
fprintf(fp, " %d-%d-%d-%d (%g)",
idihs[i].ai()+1, idihs[i].aj()+1,
}
}
-static void print_param(FILE *fp, int ftype, int i, t_param *param)
+static void printInteractionType(FILE *fp, int ftype, int i, const InteractionType &type)
{
static int pass = 0;
static int prev_ftype = NOTSET;
static int prev_i = NOTSET;
- int j;
if ( (ftype != prev_ftype) || (i != prev_i) )
{
}
fprintf(fp, "(%d) plist[%s].param[%d]",
pass, interaction_function[ftype].name, i);
- for (j = 0; j < NRFP(ftype); j++)
+ gmx::ArrayRef<const real> 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++;
}
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 )
{
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)
{
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)
{
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)
{
/ ( 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)
{
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;
static bool
-calc_vsite4fdn_param(t_param *param,
+calc_vsite4fdn_param(InteractionType *param,
int nrbond, t_mybonded *bonds,
int nrang, t_mybonded *angles)
{
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;
int set_vsites(bool bVerbose, t_atoms *atoms, PreprocessingAtomTypes *atypes,
gmx::ArrayRef<InteractionTypeParameters> 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;
{
if (interaction_function[ftype].flags & IF_VSITE)
{
- nvsite += plist[ftype].nr;
+ nvsite += plist[ftype].size();
if (ftype == F_VSITEN)
{
}
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<const real> 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)
{
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 */
{
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)
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 */
}
static void check_vsite_constraints(gmx::ArrayRef<InteractionTypeParameters> 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<const int> 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++;
}
}
static void clean_vsite_bonds(gmx::ArrayRef<InteractionTypeParameters> 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)
}
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<const int> 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++;
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
{
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;
}
}
}
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;
}
}
}
* 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;
}
}
}
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)
{
nOut, interaction_function[cftype].longname,
interaction_function[F_VSITE3OUT].longname );
}
- ps->nr = kept_i;
}
static void clean_vsite_angles(gmx::ArrayRef<InteractionTypeParameters> 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<const int> 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++;
{
/* 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
{
/* 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;
}
}
}
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;
}
}
}
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<InteractionTypeParameters> 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<const int> 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
{
/* 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;
}
}
}
/* 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<InteractionTypeParameters> plist, int natoms, bool bRmVSiteBds)
{
- int i, k, nvsite, ftype, vsite, parnr;
+ int nvsite, vsite;
int *vsite_type;
t_pindex *pindex;
at2vsitecon_t *at2vc;
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;
}
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++;
}
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
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;
}
}
/* 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 ) )
}
}
/* 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)
{
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<real, MAXFORCEPARAM> forceParam = {0.0};
if (bPBC)
{
set_pbc(&pbc, -1, box);
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<int> atoms = {i, j};
+ add_param_to_list (bond, InteractionType(atoms, forceParam));
nbond[i]++;
nbond[j]++;
}
double cc;
char buf[32];
- for (int i = 0; (i < plist->nr); i++)
+ for (auto ¶m : plist->interactionTypes)
{
if (!bParam)
{
{
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)
{
}
}
GMX_ASSERT(nrfp <= MAXFORCEPARAM/2, "Only 6 parameters may be used for an interaction");
+ std::array<real, MAXFORCEPARAM> 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);
}
}
{
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);
}
}
static void print_pl(FILE *fp, gmx::ArrayRef<const InteractionTypeParameters> 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<const int> atoms = param.atoms();
+ gmx::ArrayRef<const real> 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");
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);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff