#include "gromacs/fileio/warninp.h"
#include "gromacs/gmxpreprocess/gpp_atomtype.h"
#include "gromacs/gmxpreprocess/gpp_bond_atomtype.h"
+#include "gromacs/gmxpreprocess/grompp_impl.h"
#include "gromacs/gmxpreprocess/notset.h"
#include "gromacs/gmxpreprocess/readir.h"
#include "gromacs/gmxpreprocess/topdirs.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/stringutil.h"
-void generate_nbparams(int comb, int ftype, t_params *plist, gpp_atomtype_t atype,
- warninp_t wi)
+void generate_nbparams(int comb,
+ int ftype,
+ InteractionsOfType *interactions,
+ PreprocessingAtomTypes *atypes,
+ warninp *wi)
{
- int i, j, k = -1, nf;
int nr, nrfp;
real c, bi, bj, ci, cj, ci0, ci1, ci2, cj0, cj1, cj2;
/* Lean mean shortcuts */
- nr = get_atomtype_ntypes(atype);
+ nr = atypes->size();
nrfp = NRFP(ftype);
- snew(plist->param, nr*nr);
- plist->nr = nr*nr;
+ interactions->interactionTypes.clear();
+ std::array<real, MAXFORCEPARAM> forceParam = {NOTSET};
/* Fill the matrix with force parameters */
switch (ftype)
{
{
case eCOMB_GEOMETRIC:
/* Gromos rules */
- for (i = k = 0; (i < nr); i++)
+ for (int i = 0; (i < nr); i++)
{
- for (j = 0; (j < nr); j++, k++)
+ for (int j = 0; (j < nr); j++)
{
- for (nf = 0; (nf < nrfp); nf++)
+ for (int nf = 0; (nf < nrfp); nf++)
{
- ci = get_atomtype_nbparam(i, nf, atype);
- cj = get_atomtype_nbparam(j, nf, atype);
- 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;
}
+ interactions->interactionTypes.emplace_back(InteractionOfType({}, forceParam));
}
}
break;
case eCOMB_ARITHMETIC:
/* c0 and c1 are sigma and epsilon */
- for (i = k = 0; (i < nr); i++)
+ for (int i = 0; (i < nr); i++)
{
- for (j = 0; (j < nr); j++, k++)
+ for (int j = 0; (j < nr); j++)
{
- ci0 = get_atomtype_nbparam(i, 0, atype);
- cj0 = get_atomtype_nbparam(j, 0, atype);
- ci1 = get_atomtype_nbparam(i, 1, atype);
- cj1 = get_atomtype_nbparam(j, 1, atype);
- plist->param[k].c[0] = (fabs(ci0) + fabs(cj0))*0.5;
+ ci0 = atypes->atomNonBondedParamFromAtomType(i, 0);
+ cj0 = atypes->atomNonBondedParamFromAtomType(j, 0);
+ ci1 = atypes->atomNonBondedParamFromAtomType(i, 1);
+ cj1 = atypes->atomNonBondedParamFromAtomType(j, 1);
+ 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);
+ interactions->interactionTypes.emplace_back(InteractionOfType({}, forceParam));
}
}
break;
case eCOMB_GEOM_SIG_EPS:
/* c0 and c1 are sigma and epsilon */
- for (i = k = 0; (i < nr); i++)
+ for (int i = 0; (i < nr); i++)
{
- for (j = 0; (j < nr); j++, k++)
+ for (int j = 0; (j < nr); j++)
{
- ci0 = get_atomtype_nbparam(i, 0, atype);
- cj0 = get_atomtype_nbparam(j, 0, atype);
- ci1 = get_atomtype_nbparam(i, 1, atype);
- cj1 = get_atomtype_nbparam(j, 1, atype);
- plist->param[k].c[0] = std::sqrt(std::fabs(ci0*cj0));
+ ci0 = atypes->atomNonBondedParamFromAtomType(i, 0);
+ cj0 = atypes->atomNonBondedParamFromAtomType(j, 0);
+ ci1 = atypes->atomNonBondedParamFromAtomType(i, 1);
+ cj1 = atypes->atomNonBondedParamFromAtomType(j, 1);
+ 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);
+ interactions->interactionTypes.emplace_back(InteractionOfType({}, 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 (i = k = 0; (i < nr); i++)
+ for (int i = 0; (i < nr); i++)
{
- for (j = 0; (j < nr); j++, k++)
+ for (int j = 0; (j < nr); j++)
{
- ci0 = get_atomtype_nbparam(i, 0, atype);
- cj0 = get_atomtype_nbparam(j, 0, atype);
- ci2 = get_atomtype_nbparam(i, 2, atype);
- cj2 = get_atomtype_nbparam(j, 2, atype);
- bi = get_atomtype_nbparam(i, 1, atype);
- bj = get_atomtype_nbparam(j, 1, atype);
- plist->param[k].c[0] = std::sqrt(ci0 * cj0);
+ ci0 = atypes->atomNonBondedParamFromAtomType(i, 0);
+ cj0 = atypes->atomNonBondedParamFromAtomType(j, 0);
+ ci2 = atypes->atomNonBondedParamFromAtomType(i, 2);
+ cj2 = atypes->atomNonBondedParamFromAtomType(j, 2);
+ bi = atypes->atomNonBondedParamFromAtomType(i, 1);
+ bj = atypes->atomNonBondedParamFromAtomType(j, 1);
+ 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);
+ interactions->interactionTypes.emplace_back(InteractionOfType({}, forceParam));
}
}
}
}
-static void realloc_nb_params(gpp_atomtype_t at,
+/*! \brief Used to temporarily store the explicit non-bonded parameter
+ * combinations, which will be copied to InteractionsOfType. */
+struct t_nbparam
+{
+ //! Has this combination been set.
+ bool bSet;
+ //! The non-bonded parameters
+ real c[4];
+};
+
+static void realloc_nb_params(PreprocessingAtomTypes *atypes,
t_nbparam ***nbparam, t_nbparam ***pair)
{
/* Add space in the non-bonded parameters matrix */
- int atnr = get_atomtype_ntypes(at);
+ int atnr = atypes->size();
srenew(*nbparam, atnr);
snew((*nbparam)[atnr-1], atnr);
if (pair)
}
}
+int copy_nbparams(t_nbparam **param, int ftype, InteractionsOfType *interactions, int nr)
+{
+ int nrfp, ncopy;
+
+ nrfp = NRFP(ftype);
+
+ ncopy = 0;
+ for (int i = 0; i < nr; i++)
+ {
+ for (int j = 0; j <= i; j++)
+ {
+ GMX_RELEASE_ASSERT(param, "Must have valid parameters");
+ if (param[i][j].bSet)
+ {
+ for (int f = 0; f < nrfp; f++)
+ {
+ interactions->interactionTypes[nr*i+j].setForceParameter(f, param[i][j].c[f]);
+ interactions->interactionTypes[nr*j+i].setForceParameter(f, param[i][j].c[f]);
+ }
+ ncopy++;
+ }
+ }
+ }
+
+ return ncopy;
+}
+
+void free_nbparam(t_nbparam **param, int nr)
+{
+ int i;
+
+ GMX_RELEASE_ASSERT(param, "Must have valid parameters");
+ for (i = 0; i < nr; i++)
+ {
+ GMX_RELEASE_ASSERT(param[i], "Must have valid parameters");
+ sfree(param[i]);
+ }
+ sfree(param);
+}
+
static void copy_B_from_A(int ftype, double *c)
{
int nrfpA, nrfpB, i;
}
}
-void push_at (t_symtab *symtab, gpp_atomtype_t at, t_bond_atomtype bat,
+void push_at (t_symtab *symtab, PreprocessingAtomTypes *at, PreprocessingBondAtomType *bondAtomType,
char *line, int nb_funct,
t_nbparam ***nbparam, t_nbparam ***pair,
- warninp_t wi)
+ warninp *wi)
{
typedef struct {
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",
auto message = gmx::formatString("Invalid function type %d in push_at", nb_funct);
warning_error_and_exit(wi, message, FARGS);
}
- for (j = nfp0; (j < MAXFORCEPARAM); j++)
+ for (int j = nfp0; (j < MAXFORCEPARAM); j++)
{
c[j] = 0.0;
}
+ 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];
}
- if ((batype_nr = get_bond_atomtype_type(btype, bat)) == NOTSET)
- {
- add_bond_atomtype(bat, symtab, btype);
- }
- batype_nr = get_bond_atomtype_type(btype, bat);
+ InteractionOfType interactionType({}, forceParam, "");
+
+ batype_nr = bondAtomType->addBondAtomType(symtab, btype);
- if ((nr = get_atomtype_type(type, at)) != NOTSET)
+ if ((nr = at->atomTypeFromName(type)) != NOTSET)
{
auto message = gmx::formatString
("Atomtype %s was defined previously (e.g. in the forcefield files), "
"you should override the previous definition, then you could choose "
"to suppress this warning with -maxwarn.", type);
warning(wi, message);
- if ((nr = set_atomtype(nr, at, symtab, atom, type, param, batype_nr,
- atomnr)) == NOTSET)
+ if ((nr = at->setType(nr, symtab, *atom, type, interactionType, batype_nr,
+ atomnr)) == NOTSET)
{
auto message = gmx::formatString("Replacing atomtype %s failed", type);
warning_error_and_exit(wi, message, FARGS);
}
}
- else if ((add_atomtype(at, symtab, atom, type, param,
- batype_nr, atomnr)) == NOTSET)
+ else if ((at->addType(symtab, *atom, type, interactionType,
+ batype_nr, atomnr)) == NOTSET)
{
auto message = gmx::formatString("Adding atomtype %s failed", type);
warning_error_and_exit(wi, message, FARGS);
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(t_params * bt,
- const t_param * b,
- int nral,
- int ftype,
- bool bAllowRepeat,
- const char * line,
- warninp_t wi)
+static void push_bondtype(InteractionsOfType *bt,
+ const InteractionOfType &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(
+ InteractionOfType(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->interactionTypes.emplace_back(InteractionOfType(atoms, forceParam, b.interactionTypeName()));
+ }
+}
+
+static std::vector<int> atomTypesFromAtomNames(const PreprocessingAtomTypes *atomTypes,
+ const PreprocessingBondAtomType *bondAtomTypes,
+ gmx::ArrayRef<const char[20]> atomNames,
+ warninp *wi)
+{
+
+ GMX_RELEASE_ASSERT(!(!atomNames.empty() && !atomTypes && !bondAtomTypes),
+ "Need to have either valid atomtypes or bondatomtypes object");
- bt->nr += 2;
+ std::vector<int> atomTypesFromAtomNames;
+ for (const auto &name : atomNames)
+ {
+ if (atomTypes != nullptr)
+ {
+ int atomType = atomTypes->atomTypeFromName(name);
+ if (atomType == NOTSET)
+ {
+ auto message = gmx::formatString("Unknown atomtype %s\n", name);
+ warning_error_and_exit(wi, message, FARGS);
+ }
+ atomTypesFromAtomNames.emplace_back(atomType);
+ }
+ else if (bondAtomTypes != nullptr)
+ {
+ int atomType = bondAtomTypes->bondAtomTypeFromName(name);
+ if (atomType == NOTSET)
+ {
+ auto message = gmx::formatString("Unknown bond_atomtype %s\n", name);
+ warning_error_and_exit(wi, message, FARGS);
+ }
+ atomTypesFromAtomNames.emplace_back(atomType);
+ }
}
+ return atomTypesFromAtomNames;
}
-void push_bt(directive d, t_params bt[], int nral,
- gpp_atomtype_t at,
- t_bond_atomtype bat, char *line,
- warninp_t wi)
+
+void push_bt(Directive d,
+ gmx::ArrayRef<InteractionsOfType> bt,
+ int nral,
+ PreprocessingAtomTypes *at,
+ PreprocessingBondAtomType *bondAtomType,
+ 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))
+ if ((bondAtomType && at) || (!bondAtomType && !at))
{
- gmx_incons("You should pass either bat or at to push_bt");
+ gmx_incons("You should pass either bondAtomType or at to push_bt");
}
/* Make format string (nral ints+functype) */
}
}
}
- for (i = 0; (i < nral); i++)
- {
- if (at && ((p.a[i] = get_atomtype_type(alc[i], at)) == 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))
- {
- auto message = gmx::formatString("Unknown bond_atomtype %s\n", alc[i]);
- warning_error_and_exit(wi, message, FARGS);
- }
- }
- for (i = 0; (i < nrfp); i++)
+ std::vector<int> atomTypes = atomTypesFromAtomNames(at,
+ bondAtomType,
+ gmx::arrayRefFromArray(alc, nral),
+ wi);
+ std::array<real, MAXFORCEPARAM> forceParam;
+ 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]), InteractionOfType(atomTypes, forceParam), nral, ftype, FALSE, line, wi);
}
-void push_dihedraltype(directive d, t_params bt[],
- t_bond_atomtype bat, char *line,
- warninp_t wi)
+void push_dihedraltype(Directive d, gmx::ArrayRef<InteractionsOfType> bt,
+ PreprocessingBondAtomType *bondAtomType, 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 = bondAtomType->bondAtomTypeFromName(alc[i])) == 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]), InteractionOfType(atoms, forceParam), 4, ftype, bAllowRepeat, line, wi);
}
-void push_nbt(directive d, t_nbparam **nbt, gpp_atomtype_t atype,
+void push_nbt(Directive d, t_nbparam **nbt, PreprocessingAtomTypes *atypes,
char *pline, int nb_funct,
- warninp_t wi)
+ warninp *wi)
{
/* swap the atoms */
const char *form3 = "%*s%*s%*s%lf%lf%lf";
}
/* Put the parameters in the matrix */
- if ((ai = get_atomtype_type (a0, atype)) == NOTSET)
+ if ((ai = atypes->atomTypeFromName(a0)) == NOTSET)
{
auto message = gmx::formatString("Atomtype %s not found", a0);
warning_error_and_exit(wi, message, FARGS);
}
- if ((aj = get_atomtype_type (a1, atype)) == NOTSET)
+ if ((aj = atypes->atomTypeFromName(a1)) == NOTSET)
{
auto message = gmx::formatString("Atomtype %s not found", a1);
warning_error_and_exit(wi, message, FARGS);
}
void
-push_cmaptype(directive d, t_params bt[], int nral, gpp_atomtype_t at,
- t_bond_atomtype bat, char *line,
- warninp_t wi)
+push_cmaptype(Directive d,
+ gmx::ArrayRef<InteractionsOfType> bt,
+ int nral,
+ PreprocessingAtomTypes *atomtypes,
+ PreprocessingBondAtomType *bondAtomType,
+ 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 i, 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;
nrfpB = strtol(alc[7], nullptr, 10)*strtol(alc[7], nullptr, 10);
nrfp = nrfpA+nrfpB;
- /* Allocate memory for the CMAP grid */
- bt[F_CMAP].ncmap += nrfp;
- srenew(bt[F_CMAP].cmap, bt[F_CMAP].ncmap);
-
/* Read in CMAP parameters */
sl = 0;
- for (i = 0; i < ncmap; i++)
+ for (int i = 0; i < ncmap; i++)
{
while (isspace(*(line+start+sl)))
{
}
nn = sscanf(line+start+sl, " %s ", s);
sl += strlen(s);
- bt[F_CMAP].cmap[i+(bt[F_CMAP].ncmap)-nrfp] = strtod(s, nullptr);
+ bt[F_CMAP].cmap.emplace_back(strtod(s, nullptr));
if (nn == 1)
{
/* Check do that we got the number of parameters we expected */
if (read_cmap == nrfpA)
{
- for (i = 0; i < ncmap; i++)
+ for (int i = 0; i < ncmap; i++)
{
- bt[F_CMAP].cmap[i+ncmap] = bt[F_CMAP].cmap[i];
+ bt[F_CMAP].cmap.emplace_back(bt[F_CMAP].cmap[i]);
}
}
else
/* Set grid spacing and the number of grids (we assume these numbers to be the same for all grids
* so we can safely assign them each time
*/
- bt[F_CMAP].grid_spacing = nxcmap; /* Or nycmap, they need to be equal */
- bt[F_CMAP].nc = bt[F_CMAP].nc + 1; /* Since we are incrementing here, we need to subtract later, see (*****) */
- nct = (nral+1) * bt[F_CMAP].nc;
+ bt[F_CMAP].cmakeGridSpacing = nxcmap; /* Or nycmap, they need to be equal */
+ bt[F_CMAP].cmapAngles++; /* Since we are incrementing here, we need to subtract later, see (*****) */
+ nct = (nral+1) * bt[F_CMAP].cmapAngles;
- /* Allocate memory for the cmap_types information */
- srenew(bt[F_CMAP].cmap_types, nct);
-
- for (i = 0; (i < nral); i++)
+ for (int i = 0; (i < nral); i++)
{
- if (at && ((p.a[i] = get_bond_atomtype_type(alc[i], bat)) == NOTSET))
- {
- 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))
- {
- auto message = gmx::formatString("Unknown bond_atomtype %s\n", alc[i]);
- warning_error(wi, message);
- }
-
/* Assign a grid number to each cmap_type */
- bt[F_CMAP].cmap_types[bt[F_CMAP].nct++] = get_bond_atomtype_type(alc[i], bat);
+ GMX_RELEASE_ASSERT(bondAtomType != nullptr, "Need valid PreprocessingBondAtomType object");
+ bt[F_CMAP].cmapAtomTypes.emplace_back(bondAtomType->bondAtomTypeFromName(alc[i]));
}
/* Assign a type number to this cmap */
- bt[F_CMAP].cmap_types[bt[F_CMAP].nct++] = bt[F_CMAP].nc-1; /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */
+ bt[F_CMAP].cmapAtomTypes.emplace_back(bt[F_CMAP].cmapAngles-1); /* Since we inremented earlier, we need to subtrac here, to get the types right (****) */
/* Check for the correct number of atoms (again) */
- if (bt[F_CMAP].nct != nct)
+ if (bt[F_CMAP].nct() != nct)
{
- auto message = gmx::formatString("Incorrect number of atom types (%d) in cmap type %d\n", nct, bt[F_CMAP].nc);
+ auto message = gmx::formatString("Incorrect number of atom types (%d) in cmap type %d\n", nct, bt[F_CMAP].cmapAngles);
warning_error(wi, message);
}
-
- /* Is this correct?? */
- for (i = 0; i < MAXFORCEPARAM; i++)
- {
- p.c[i] = NOTSET;
- }
+ std::vector<int> atomTypes = atomTypesFromAtomNames(atomtypes,
+ bondAtomType,
+ gmx::constArrayRefFromArray(alc, nral),
+ wi);
+ 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]), InteractionOfType(atomTypes, forceParam), nral, ftype, FALSE, line, wi);
}
char *resnumberic,
char *resname, char *name, real m0, real q0,
int typeB, char *ctypeB, real mB, real qB,
- warninp_t wi)
+ warninp *wi)
{
int j, resind = 0, resnr;
unsigned char ric;
}
void push_atom(t_symtab *symtab, t_block *cgs,
- t_atoms *at, gpp_atomtype_t atype, char *line, int *lastcg,
- warninp_t wi)
+ t_atoms *at, PreprocessingAtomTypes *atypes, char *line, int *lastcg,
+ warninp *wi)
{
int nr, ptype;
int cgnumber, atomnr, type, typeB, nscan;
return;
}
sscanf(id, "%d", &atomnr);
- if ((type = get_atomtype_type(ctype, atype)) == NOTSET)
+ if ((type = atypes->atomTypeFromName(ctype)) == NOTSET)
{
auto message = gmx::formatString("Atomtype %s not found", ctype);
warning_error_and_exit(wi, message, FARGS);
}
- ptype = get_atomtype_ptype(type, atype);
+ ptype = atypes->atomParticleTypeFromAtomType(type);
/* Set default from type */
- q0 = get_atomtype_qA(type, atype);
- m0 = get_atomtype_massA(type, atype);
+ q0 = atypes->atomChargeFromAtomType(type);
+ m0 = atypes->atomMassFromAtomType(type);
typeB = type;
qB = q0;
mB = m0;
m0 = mB = m;
if (nscan > 2)
{
- if ((typeB = get_atomtype_type(ctypeB, atype)) == NOTSET)
+ if ((typeB = atypes->atomTypeFromName(ctypeB)) == NOTSET)
{
auto message = gmx::formatString("Atomtype %s not found", ctypeB);
warning_error_and_exit(wi, message, FARGS);
}
- qB = get_atomtype_qA(typeB, atype);
- mB = get_atomtype_massA(typeB, atype);
+ qB = atypes->atomChargeFromAtomType(typeB);
+ mB = atypes->atomMassFromAtomType(typeB);
if (nscan > 3)
{
qB = qb;
push_cg(cgs, lastcg, cgnumber, nr);
- push_atom_now(symtab, at, atomnr, get_atomtype_atomnumber(type, atype),
+ push_atom_now(symtab, at, atomnr, atypes->atomNumberFromAtomType(type),
type, ctype, ptype, resnumberic,
resname, name, m0, q0, typeB,
typeB == type ? ctype : ctypeB, mB, qB, wi);
}
-void push_molt(t_symtab *symtab, int *nmol, t_molinfo **mol, char *line,
- warninp_t wi)
+void push_molt(t_symtab *symtab,
+ std::vector<MoleculeInformation> *mol,
+ char *line,
+ warninp *wi)
{
char type[STRLEN];
- int nrexcl, i;
- t_molinfo *newmol;
+ int nrexcl;
if ((sscanf(line, "%s%d", type, &nrexcl)) != 2)
{
}
/* Test if this moleculetype overwrites another */
- i = 0;
- while (i < *nmol)
+ const auto found = std::find_if(mol->begin(), mol->end(),
+ [&type](const auto &m)
+ { return strcmp(*(m.name), type) == 0; });
+ if (found != mol->end())
{
- if (strcmp(*((*mol)[i].name), type) == 0)
- {
- auto message = gmx::formatString("moleculetype %s is redefined", type);
- warning_error_and_exit(wi, message, FARGS);
- }
- i++;
+ auto message = gmx::formatString("moleculetype %s is redefined", type);
+ warning_error_and_exit(wi, message, FARGS);
}
- (*nmol)++;
- srenew(*mol, *nmol);
- newmol = &((*mol)[*nmol-1]);
- init_molinfo(newmol);
+ mol->emplace_back();
+ mol->back().initMolInfo();
/* Fill in the values */
- newmol->name = put_symtab(symtab, type);
- newmol->nrexcl = nrexcl;
- newmol->excl_set = FALSE;
+ mol->back().name = put_symtab(symtab, type);
+ mol->back().nrexcl = nrexcl;
+ mol->back().excl_set = false;
}
-static bool default_nb_params(int ftype, t_params bt[], t_atoms *at,
- t_param *p, int c_start, bool bB, bool bGenPairs)
+static bool findIfAllNBAtomsMatch(gmx::ArrayRef<const int> atomsFromParameterArray,
+ gmx::ArrayRef<const int> atomsFromCurrentParameter,
+ const t_atoms *at,
+ bool bB)
{
- 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;
+ 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<InteractionsOfType> bt, t_atoms *at,
+ InteractionOfType *p, int c_start, bool bB, bool bGenPairs)
+{
+ int ti, tj, ntype;
+ bool bFound;
+ InteractionOfType *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(t_params bondtype[],
- t_atoms *at, gpp_atomtype_t atype,
- t_param *p, bool bB,
+static bool default_cmap_params(gmx::ArrayRef<InteractionsOfType> bondtype,
+ t_atoms *at, PreprocessingAtomTypes *atypes,
+ InteractionOfType *p, bool bB,
int *cmap_type, int *nparam_def,
- warninp_t wi)
+ warninp *wi)
{
- int i, nparam_found;
+ int nparam_found;
int ct;
- bool bFound = FALSE;
+ bool bFound = false;
nparam_found = 0;
ct = 0;
/* Match the current cmap angle against the list of cmap_types */
- for (i = 0; i < bondtype[F_CMAP].nct && !bFound; i += 6)
+ for (int i = 0; i < bondtype[F_CMAP].nct() && !bFound; i += 6)
{
if (bB)
{
else
{
if (
- (get_atomtype_batype(at->atom[p->a[0]].type, atype) == bondtype[F_CMAP].cmap_types[i]) &&
- (get_atomtype_batype(at->atom[p->a[1]].type, atype) == bondtype[F_CMAP].cmap_types[i+1]) &&
- (get_atomtype_batype(at->atom[p->a[2]].type, atype) == bondtype[F_CMAP].cmap_types[i+2]) &&
- (get_atomtype_batype(at->atom[p->a[3]].type, atype) == bondtype[F_CMAP].cmap_types[i+3]) &&
- (get_atomtype_batype(at->atom[p->a[4]].type, atype) == bondtype[F_CMAP].cmap_types[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;
- ct = bondtype[F_CMAP].cmap_types[i+5];
+ bFound = true;
+ ct = bondtype[F_CMAP].cmapAtomTypes[i+5];
nparam_found = 1;
}
}
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 InteractionOfType &pi,
int type_i, int type_j, int type_k, int type_l,
- const gpp_atomtype* atype)
+ const PreprocessingAtomTypes* atypes)
{
- if ((pi->ai() == -1 || get_atomtype_batype(type_i, atype) == pi->ai()) &&
- (pi->aj() == -1 || get_atomtype_batype(type_j, atype) == pi->aj()) &&
- (pi->ak() == -1 || get_atomtype_batype(type_k, atype) == pi->ak()) &&
- (pi->al() == -1 || get_atomtype_batype(type_l, atype) == pi->al()))
+ if ((pi.ai() == -1 || atypes->bondAtomTypeFromAtomType(type_i) == pi.ai()) &&
+ (pi.aj() == -1 || atypes->bondAtomTypeFromAtomType(type_j) == pi.aj()) &&
+ (pi.ak() == -1 || atypes->bondAtomTypeFromAtomType(type_k) == pi.ak()) &&
+ (pi.al() == -1 || atypes->bondAtomTypeFromAtomType(type_l) == pi.al()))
{
return
- (pi->ai() == -1 ? 0 : 1) +
- (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 default_params(int ftype, t_params bt[],
- t_atoms *at, gpp_atomtype_t atype,
- t_param *p, bool bB,
- t_param **param_def,
- int *nparam_def)
+static int findNumberOfDihedralAtomMatches(const InteractionOfType ¤tParamFromParameterArray,
+ const InteractionOfType ¶meterToAdd,
+ const t_atoms *at,
+ const PreprocessingAtomTypes *atypes,
+ bool bB)
+{
+ 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<InteractionOfType>::iterator
+defaultInteractionsOfType(int ftype, gmx::ArrayRef<InteractionsOfType> bt,
+ t_atoms *at, PreprocessingAtomTypes *atypes,
+ const InteractionOfType &p, bool bB,
+ int *nparam_def)
{
- 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;
+ 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)
+ 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 = natom_match(pt, at->atom[p->ai()].typeB, at->atom[p->aj()].typeB, at->atom[p->ak()].typeB, at->atom[p->al()].typeB, atype);
- }
- else
- {
- nmatch = natom_match(pt, at->atom[p->ai()].type, at->atom[p->aj()].type, at->atom[p->ak()].type, at->atom[p->al()].type, atype);
- }
- if (nmatch > nmatch_max)
- {
- 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;
+ //Advance iterator (like std::advance) without incrementing past end (UB)
+ const auto safeAdvance = [](auto &it, auto n, auto end) {
+ it = end-it > n ? it+n : end;
+ };
+ /* Continue from current iterator position */
+ auto nextPos = prevPos;
+ const auto endIter = bt[ftype].interactionTypes.end();
+ safeAdvance(nextPos, 2, endIter);
+ for (; nextPos < endIter && bSame; safeAdvance(nextPos, 2, endIter))
{
- 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) &&
- (get_atomtype_batype(at->atom[p->a[j]].typeB, atype) == pi->a[j])); j++)
- {
- ;
- }
- }
- else
- {
- for (j = 0; ((j < nral) &&
- (get_atomtype_batype(at->atom[p->a[j]].type, atype) == 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;
}
-void push_bond(directive d, t_params bondtype[], t_params bond[],
- t_atoms *at, gpp_atomtype_t atype, char *line,
+void push_bond(Directive d, gmx::ArrayRef<InteractionsOfType> bondtype,
+ gmx::ArrayRef<InteractionsOfType> bond,
+ t_atoms *at, PreprocessingAtomTypes *atypes, char *line,
bool bBonded, bool bGenPairs, real fudgeQQ,
bool bZero, bool *bWarn_copy_A_B,
- warninp_t wi)
+ 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++)
- {
- 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);
}
+ /* 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!
*/
+ InteractionOfType param(atoms, forceParam, "");
+ std::vector<InteractionOfType>::iterator foundAParameter = bondtype[ftype].interactionTypes.end();
+ std::vector<InteractionOfType>::iterator foundBParameter = bondtype[ftype].interactionTypes.end();
if (bBonded)
{
- bFoundA = default_params(ftype, bondtype, at, atype, ¶m, FALSE, ¶m_defA, &nparam_defA);
- if (bFoundA)
+ foundAParameter = defaultInteractionsOfType(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;
}
- bFoundB = default_params(ftype, bondtype, at, atype, ¶m, TRUE, ¶m_defB, &nparam_defB);
- if (bFoundB)
+ else if (NRFPA(ftype) == 0)
+ {
+ bFoundA = true;
+ }
+ foundBParameter = defaultInteractionsOfType(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);
}
}
}
-void push_cmap(directive d, t_params bondtype[], t_params bond[],
- t_atoms *at, gpp_atomtype_t atype, char *line,
- warninp_t wi)
+void push_cmap(Directive d, gmx::ArrayRef<InteractionsOfType> bondtype,
+ gmx::ArrayRef<InteractionsOfType> bond,
+ t_atoms *at, PreprocessingAtomTypes *atypes, char *line,
+ warninp *wi)
{
const char *aaformat[MAXATOMLIST+1] =
{
"%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++)
+ std::vector<int> atoms;
+ for (int j = 0; (j < nral); j++)
{
- param.a[j] = aa[j]-1;
+ atoms.emplace_back(aa[j]-1);
}
- for (j = 0; (j < MAXFORCEPARAM); j++)
- {
- param.c[j] = 0.0;
- }
-
+ std::array<real, MAXFORCEPARAM> forceParam = {0.0};
+ InteractionOfType param(atoms, forceParam, "");
/* Get the cmap type for this cmap angle */
- bFound = default_cmap_params(bondtype, at, atype, ¶m, FALSE, &cmap_type, &ncmap_params, wi);
+ bFound = default_cmap_params(bondtype, at, atypes, ¶m, FALSE, &cmap_type, &ncmap_params, wi);
/* We want exactly one parameter (the cmap type in state A (currently no state B) back */
if (bFound && ncmap_params == 1)
{
/* 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);
}
}
-void push_vsitesn(directive d, t_params bond[],
+void push_vsitesn(Directive d, gmx::ArrayRef<InteractionsOfType> bond,
t_atoms *at, char *line,
- warninp_t wi)
+ 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], InteractionOfType(atoms, forceParam));
}
sfree(atc);
sfree(weight);
}
-void push_mol(int nrmols, t_molinfo mols[], char *pline, int *whichmol,
+void push_mol(gmx::ArrayRef<MoleculeInformation> mols, char *pline, int *whichmol,
int *nrcopies,
- warninp_t wi)
+ warninp *wi)
{
char type[STRLEN];
int nrci = 0;
int matchci = -1;
int matchcs = -1;
- for (int i = 0; i < nrmols; i++)
+ int i = 0;
+ for (const auto &mol : mols)
{
- if (strcmp(type, *(mols[i].name)) == 0)
+ if (strcmp(type, *(mol.name)) == 0)
{
nrcs++;
matchcs = i;
}
- if (gmx_strcasecmp(type, *(mols[i].name)) == 0)
+ if (gmx_strcasecmp(type, *(mol.name)) == 0)
{
nrci++;
matchci = i;
}
+ i++;
}
if (nrcs == 1)
}
}
-void push_excl(char *line, gmx::ExclusionBlocks *b2, warninp_t wi)
+void push_excl(char *line, gmx::ArrayRef<gmx::ExclusionBlock> b2, warninp *wi)
{
int i, j;
int n;
return;
}
- if ((1 <= i) && (i <= b2->nr))
+ if ((1 <= i) && (i <= b2.ssize()))
{
i--;
}
n = sscanf(line, format, &j);
if (n == 1)
{
- if ((1 <= j) && (j <= b2->nr))
+ if ((1 <= j) && (j <= b2.ssize()))
{
j--;
- srenew(b2->a[i], ++(b2->nra[i]));
- b2->a[i][b2->nra[i]-1] = j;
+ b2[i].atomNumber.push_back(j);
/* also add the reverse exclusion! */
- srenew(b2->a[j], ++(b2->nra[j]));
- b2->a[j][b2->nra[j]-1] = i;
+ b2[j].atomNumber.push_back(i);
strcat(base, "%*d");
}
else
{
- auto message = gmx::formatString("Invalid Atomnr j: %d, b2->nr: %d\n", j, b2->nr);
+ auto message = gmx::formatString("Invalid Atomnr j: %d, b2->nr: %zu\n", j, b2.size());
warning_error_and_exit(wi, message, FARGS);
}
}
while (n == 1);
}
-int add_atomtype_decoupled(t_symtab *symtab, gpp_atomtype_t at,
+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 = add_atomtype(at, symtab, &atom, "decoupled", ¶m, -1, 0);
+ std::array<real, MAXFORCEPARAM> forceParam = {0.0};
+ nr = at->addType(symtab, atom, "decoupled", InteractionOfType({}, forceParam, ""), -1, 0);
/* Add space in the non-bonded parameters matrix */
realloc_nb_params(at, nbparam, pair);
return nr;
}
-static void convert_pairs_to_pairsQ(t_params *plist,
+static void convert_pairs_to_pairsQ(gmx::ArrayRef<InteractionsOfType> interactions,
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<InteractionOfType> paramnew;
- paramp1 = plist[F_LJ14].param;
- paramp2 = plist[F_LJC14_Q].param;
+ gmx::ArrayRef<const InteractionOfType> paramp1 = interactions[F_LJ14].interactionTypes;
+ gmx::ArrayRef<const InteractionOfType> paramp2 = interactions[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(InteractionOfType(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;
+ interactions[F_LJC14_Q].interactionTypes = paramnew;
/* Empty the LJ14 pairlist */
- plist[F_LJ14].nr = 0;
- plist[F_LJ14].param = nullptr;
+ interactions[F_LJ14].interactionTypes.clear();
}
-static void generate_LJCpairsNB(t_molinfo *mol, int nb_funct, t_params *nbp, warninp_t wi)
+static void generate_LJCpairsNB(MoleculeInformation *mol, int nb_funct, InteractionsOfType *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->interactions[F_LJC_PAIRS_NB], InteractionOfType(atoms, forceParam));
}
}
}
static void decouple_atoms(t_atoms *atoms, int atomtype_decouple,
int couple_lam0, int couple_lam1,
- const char *mol_name, warninp_t wi)
+ const char *mol_name, warninp *wi)
{
int i;
}
}
-void convert_moltype_couple(t_molinfo *mol, int atomtype_decouple, real fudgeQQ,
+void convert_moltype_couple(MoleculeInformation *mol, int atomtype_decouple, real fudgeQQ,
int couple_lam0, int couple_lam1,
- bool bCoupleIntra, int nb_funct, t_params *nbp,
- warninp_t wi)
+ bool bCoupleIntra, int nb_funct, InteractionsOfType *nbp,
+ warninp *wi)
{
- convert_pairs_to_pairsQ(mol->plist, fudgeQQ, &mol->atoms);
+ convert_pairs_to_pairsQ(mol->interactions, fudgeQQ, &mol->atoms);
if (!bCoupleIntra)
{
generate_LJCpairsNB(mol, nb_funct, nbp, wi);
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff