#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);
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff