static int multi_body_bondeds_count(const gmx_mtop_t *mtop)
{
- int n, nmol, ftype;
- gmx_mtop_ilistloop_t iloop;
- const t_ilist *il;
+ int n, nmol, ftype;
+ gmx_mtop_ilistloop_t iloop;
+ const InteractionList *il;
n = 0;
iloop = gmx_mtop_ilistloop_init(mtop);
if ((interaction_function[ftype].flags & IF_BOND) &&
NRAL(ftype) > 2)
{
- n += nmol*il[ftype].nr/(1 + NRAL(ftype));
+ n += nmol*il[ftype].size()/(1 + NRAL(ftype));
}
}
}
/*! \brief Walks over the constraints out from the local atoms into the non-local atoms and adds them to a list */
static void walk_out(int con, int con_offset, int a, int offset, int nrec,
- int ncon1, const t_iatom *ia1, const t_iatom *ia2,
+ gmx::ArrayRef<const int> ia1,
+ gmx::ArrayRef<const int> ia2,
const t_blocka *at2con,
const gmx_ga2la_t &ga2la, gmx_bool bHomeConnect,
gmx_domdec_constraints_t *dc,
il_local->nalloc = over_alloc_dd(il_local->nr+3);
srenew(il_local->iatoms, il_local->nalloc);
}
- iap = constr_iatomptr(ncon1, ia1, ia2, con);
+ iap = constr_iatomptr(ia1, ia2, con);
il_local->iatoms[il_local->nr++] = iap[0];
a1_gl = offset + iap[1];
a2_gl = offset + iap[2];
if (coni != con)
{
/* Walk further */
- iap = constr_iatomptr(ncon1, ia1, ia2, coni);
+ iap = constr_iatomptr(ia1, ia2, coni);
if (a == iap[1])
{
b = iap[2];
if (!ga2la.findHome(offset + b))
{
walk_out(coni, con_offset, b, offset, nrec-1,
- ncon1, ia1, ia2, at2con,
+ ia1, ia2, at2con,
ga2la, FALSE, dc, dcc, il_local, ireq);
}
}
static void atoms_to_settles(gmx_domdec_t *dd,
const gmx_mtop_t *mtop,
const int *cginfo,
- const int *const*at2settle_mt,
+ gmx::ArrayRef < const std::vector < int>> at2settle_mt,
int cg_start, int cg_end,
t_ilist *ils_local,
std::vector<int> *ireq)
if (settle >= 0)
{
- int offset = a_gl - a_mol;
+ int offset = a_gl - a_mol;
- t_iatom *ia1 = mtop->moltype[molb->type].ilist[F_SETTLE].iatoms;
+ const int *ia1 = mtop->moltype[molb->type].ilist[F_SETTLE].iatoms.data();
- int a_gls[3];
- gmx_bool bAssign = FALSE;
- int nlocal = 0;
+ int a_gls[3];
+ gmx_bool bAssign = FALSE;
+ int nlocal = 0;
for (int sa = 0; sa < nral; sa++)
{
int a_glsa = offset + ia1[settle*(1+nral)+1+sa];
std::vector<int> *ireq)
{
const t_blocka *at2con;
- int ncon1;
- t_iatom *ia1, *ia2, *iap;
int b_lo, offset, b_mol, i, con, con_offset;
gmx_domdec_constraints_t *dc = dd->constraints;
int molnr, a_mol;
mtopGetMolblockIndex(mtop, a_gl, &mb, &molnr, &a_mol);
- const gmx_molblock_t *molb = &mtop->molblock[mb];
+ const gmx_molblock_t &molb = mtop->molblock[mb];
- ncon1 = mtop->moltype[molb->type].ilist[F_CONSTR].nr/NRAL(F_SETTLE);
-
- ia1 = mtop->moltype[molb->type].ilist[F_CONSTR].iatoms;
- ia2 = mtop->moltype[molb->type].ilist[F_CONSTRNC].iatoms;
+ gmx::ArrayRef<const int> ia1 = mtop->moltype[molb.type].ilist[F_CONSTR].iatoms;
+ gmx::ArrayRef<const int> ia2 = mtop->moltype[molb.type].ilist[F_CONSTRNC].iatoms;
/* Calculate the global constraint number offset for the molecule.
* This is only required for the global index to make sure
/* The global atom number offset for this molecule */
offset = a_gl - a_mol;
- at2con = &at2con_mt[molb->type];
+ at2con = &at2con_mt[molb.type];
for (i = at2con->index[a_mol]; i < at2con->index[a_mol+1]; i++)
{
- con = at2con->a[i];
- iap = constr_iatomptr(ncon1, ia1, ia2, con);
+ con = at2con->a[i];
+ const int *iap = constr_iatomptr(ia1, ia2, con);
if (a_mol == iap[1])
{
b_mol = iap[2];
* after this first call.
*/
walk_out(con, con_offset, b_mol, offset, nrec,
- ncon1, ia1, ia2, at2con,
+ ia1, ia2, at2con,
ga2la, TRUE, dc, dcc, ilc_local, ireq);
}
}
t_ilist *ilc_local, *ils_local;
std::vector<int> *ireq;
gmx::ArrayRef<const t_blocka> at2con_mt;
- const int *const*at2settle_mt;
gmx::HashedMap<int> *ga2la_specat;
int at_end, i, j;
t_iatom *iap;
ireq = nullptr;
}
+ gmx::ArrayRef < const std::vector < int>> at2settle_mt;
+ /* When settle works inside charge groups, we assigned them already */
if (dd->bInterCGsettles)
{
// TODO Perhaps gmx_domdec_constraints_t should keep a valid constr?
at2settle_mt = constr->atom2settle_moltype();
ils_local->nr = 0;
}
- else
- {
- /* Settle works inside charge groups, we assigned them already */
- at2settle_mt = nullptr;
- }
- if (at2settle_mt == nullptr)
+ if (at2settle_mt.empty())
{
atoms_to_constraints(dd, mtop, cginfo, at2con_mt, nrec,
ilc_local, ireq);
molb = &mtop->molblock[mb];
dc->molb_con_offset[mb] = ncon;
dc->molb_ncon_mol[mb] =
- mtop->moltype[molb->type].ilist[F_CONSTR].nr/3 +
- mtop->moltype[molb->type].ilist[F_CONSTRNC].nr/3;
+ mtop->moltype[molb->type].ilist[F_CONSTR].size()/3 +
+ mtop->moltype[molb->type].ilist[F_CONSTRNC].size()/3;
ncon += molb->nmol*dc->molb_ncon_mol[mb];
}
}
/*! \brief Run the reverse ilist generation and store it in r_il when \p bAssign = TRUE */
-static int low_make_reverse_ilist(const t_ilist *il_mt, const t_atom *atom,
+static int low_make_reverse_ilist(const InteractionList *il_mt,
+ const t_atom *atom,
gmx::ArrayRef < const std::vector < int>> vsitePbc,
int *count,
gmx_bool bConstr, gmx_bool bSettle,
gmx_bool bLinkToAllAtoms,
gmx_bool bAssign)
{
- int ftype, nral, i, j, nlink, link;
- const t_ilist *il;
- const t_iatom *ia;
+ int ftype, j, nlink, link;
int a;
int nint;
- gmx_bool bVSite;
nint = 0;
for (ftype = 0; ftype < F_NRE; ftype++)
(bConstr && (ftype == F_CONSTR || ftype == F_CONSTRNC)) ||
(bSettle && ftype == F_SETTLE))
{
- bVSite = ((interaction_function[ftype].flags & IF_VSITE) != 0u);
- nral = NRAL(ftype);
- il = &il_mt[ftype];
- for (i = 0; i < il->nr; i += 1+nral)
+ const bool bVSite = ((interaction_function[ftype].flags & IF_VSITE) != 0u);
+ const int nral = NRAL(ftype);
+ const auto &il = il_mt[ftype];
+ for (int i = 0; i < il.size(); i += 1+nral)
{
- ia = il->iatoms + i;
+ const int* ia = il.iatoms.data() + i;
if (bLinkToAllAtoms)
{
if (bVSite)
}
/*! \brief Make the reverse ilist: a list of bonded interactions linked to atoms */
-static int make_reverse_ilist(const t_ilist *ilist,
+static int make_reverse_ilist(const InteractionList *ilist,
const t_atoms *atoms,
gmx::ArrayRef < const std::vector < int>> vsitePbc,
gmx_bool bConstr, gmx_bool bSettle,
atoms_global.nr = mtop->natoms;
atoms_global.atom = nullptr; /* Only used with virtual sites */
+ GMX_RELEASE_ASSERT(mtop->intermolecular_ilist.get(), "We should have an ilist when intermolecular interactions are on");
+
*nint +=
- make_reverse_ilist(mtop->intermolecular_ilist, &atoms_global,
+ make_reverse_ilist(mtop->intermolecular_ilist->data(),
+ &atoms_global,
gmx::EmptyArrayRef(),
rt.bConstr, rt.bSettle, rt.bBCheck, FALSE,
&rt.ril_intermol);
atoms.nr = mtop->natoms;
atoms.atom = nullptr;
- make_reverse_ilist(mtop->intermolecular_ilist, &atoms,
+ GMX_RELEASE_ASSERT(mtop->intermolecular_ilist.get(), "We should have an ilist when intermolecular interactions are on");
+
+ make_reverse_ilist(mtop->intermolecular_ilist->data(),
+ &atoms,
gmx::EmptyArrayRef(),
FALSE, FALSE, FALSE, TRUE, &ril_intermol);
}
{
if (dd_check_ftype(ftype, bBCheck, FALSE, FALSE))
{
- const t_ilist *il = &molt->ilist[ftype];
+ const auto *il = &molt->ilist[ftype];
int nral = NRAL(ftype);
if (nral > 1)
{
- for (int i = 0; i < il->nr; i += 1+nral)
+ for (int i = 0; i < il->size(); i += 1+nral)
{
for (int ai = 0; ai < nral; ai++)
{
}
/*! \brief Set the distance, function type and atom indices for the longest atom distance involved in intermolecular interactions for two-body and multi-body bonded interactions */
-static void bonded_distance_intermol(const t_ilist *ilists_intermol,
+static void bonded_distance_intermol(const InteractionList *ilists_intermol,
gmx_bool bBCheck,
const rvec *x, int ePBC, const matrix box,
bonded_distance_t *bd_2b,
{
if (dd_check_ftype(ftype, bBCheck, FALSE, FALSE))
{
- const t_ilist *il = &ilists_intermol[ftype];
+ const auto *il = &ilists_intermol[ftype];
int nral = NRAL(ftype);
/* No nral>1 check here, since intermol interactions always
* have nral>=2 (and the code is also correct for nral=1).
*/
- for (int i = 0; i < il->nr; i += 1+nral)
+ for (int i = 0; i < il->size(); i += 1+nral)
{
for (int ai = 0; ai < nral; ai++)
{
for (int i = 0; i < F_NRE; i++)
{
if ((interaction_function[i].flags & IF_VSITE) &&
- molt.ilist[i].nr > 0)
+ molt.ilist[i].size() > 0)
{
hasVsite = true;
}
if (moltypeHasVsite(*molt))
{
+ /* Convert to old, deprecated format */
+ t_ilist ilist[F_NRE];
+ for (int ftype = 0; ftype < F_NRE; ftype++)
+ {
+ if (interaction_function[ftype].flags & IF_VSITE)
+ {
+ ilist[ftype].nr = molt->ilist[ftype].size();
+ ilist[ftype].iatoms = const_cast<int *>(molt->ilist[ftype].iatoms.data());
+ }
+ }
+
construct_vsites(nullptr, xs, 0.0, nullptr,
- ffparams->iparams, molt->ilist,
+ ffparams->iparams, ilist,
epbcNONE, TRUE, nullptr, nullptr);
}
{
if (ir->ePBC != epbcNONE)
{
- mk_graph_ilist(nullptr, molt.ilist, 0, molt.atoms.nr, FALSE, FALSE,
- &graph);
+ mk_graph_moltype(molt, &graph);
}
std::vector<int> at2cg = make_at2cg(molt.cgs);
gmx_fatal(FARGS, "The combination of intermolecular interactions, charge groups and domain decomposition is not supported. Use cutoff-scheme=Verlet (which removes the charge groups) or run without domain decomposition.");
}
- bonded_distance_intermol(mtop->intermolecular_ilist,
+ GMX_RELEASE_ASSERT(mtop->intermolecular_ilist.get(), "We should have an ilist when intermolecular interactions are on");
+
+ bonded_distance_intermol(mtop->intermolecular_ilist->data(),
bBCheck,
x, ir->ePBC, box,
&bd_2b, &bd_mb);
int j;
std::vector<real> atomCharges;
std::vector<real> atomMasses;
- const t_ilist *ilist;
tng_bond_t tngBond;
char datatype;
{
if (IS_CHEMBOND(i))
{
- ilist = &molType->ilist[i];
- if (ilist)
+ const InteractionList &ilist = molType->ilist[i];
+ j = 1;
+ while (j < ilist.size())
{
- j = 1;
- while (j < ilist->nr)
- {
- tng_molecule_bond_add(tng, tngMol, ilist->iatoms[j], ilist->iatoms[j+1], &tngBond);
- j += 3;
- }
+ tng_molecule_bond_add(tng, tngMol, ilist.iatoms[j], ilist.iatoms[j + 1], &tngBond);
+ j += 3;
}
}
}
/* Settle is described using three atoms */
- ilist = &molType->ilist[F_SETTLE];
- if (ilist)
+ const InteractionList &ilist = molType->ilist[F_SETTLE];
+ j = 1;
+ while (j < ilist.size())
{
- j = 1;
- while (j < ilist->nr)
- {
- tng_molecule_bond_add(tng, tngMol, ilist->iatoms[j], ilist->iatoms[j+1], &tngBond);
- tng_molecule_bond_add(tng, tngMol, ilist->iatoms[j], ilist->iatoms[j+2], &tngBond);
- j += 4;
- }
+ tng_molecule_bond_add(tng, tngMol, ilist.iatoms[j], ilist.iatoms[j + 1], &tngBond);
+ tng_molecule_bond_add(tng, tngMol, ilist.iatoms[j], ilist.iatoms[j + 2], &tngBond);
+ j += 4;
}
/* First copy atom charges and masses, first atom by atom and then copy the memory for the molecule instances.
* FIXME: Atom B state data should also be written to TNG (v 2.0?) */
const t_atoms *atoms;
const t_atom *at;
const t_resinfo *resInfo;
- const t_ilist *ilist;
int nameIndex;
int atom_offset = 0;
tng_molecule_t mol, iterMol;
{
if (IS_CHEMBOND(k))
{
- ilist = &molType.ilist[k];
- if (ilist)
+ const InteractionList &ilist = molType.ilist[k];
+ for (int l = 1; l < ilist.size(); l += 3)
{
- int l = 1;
- while (l < ilist->nr)
+ int atom1, atom2;
+ atom1 = ilist.iatoms[l] + atom_offset;
+ atom2 = ilist.iatoms[l + 1] + atom_offset;
+ if (getGroupType(&mtop->groups, egcCompressedX, atom1) == 0 &&
+ getGroupType(&mtop->groups, egcCompressedX, atom2) == 0)
{
- int atom1, atom2;
- atom1 = ilist->iatoms[l] + atom_offset;
- atom2 = ilist->iatoms[l+1] + atom_offset;
- if (getGroupType(&mtop->groups, egcCompressedX, atom1) == 0 &&
- getGroupType(&mtop->groups, egcCompressedX, atom2) == 0)
- {
- tng_molecule_bond_add(tng, mol, ilist->iatoms[l],
- ilist->iatoms[l+1], &tngBond);
- }
- l += 3;
+ tng_molecule_bond_add(tng, mol, ilist.iatoms[l],
+ ilist.iatoms[l + 1], &tngBond);
}
}
}
}
/* Settle is described using three atoms */
- ilist = &molType.ilist[F_SETTLE];
- if (ilist)
+ const InteractionList &ilist = molType.ilist[F_SETTLE];
+ for (int l = 1; l < ilist.size(); l += 4)
{
- int l = 1;
- while (l < ilist->nr)
+ int atom1, atom2, atom3;
+ atom1 = ilist.iatoms[l] + atom_offset;
+ atom2 = ilist.iatoms[l + 1] + atom_offset;
+ atom3 = ilist.iatoms[l + 2] + atom_offset;
+ if (getGroupType(&mtop->groups, egcCompressedX, atom1) == 0)
{
- int atom1, atom2, atom3;
- atom1 = ilist->iatoms[l] + atom_offset;
- atom2 = ilist->iatoms[l+1] + atom_offset;
- atom3 = ilist->iatoms[l+2] + atom_offset;
- if (getGroupType(&mtop->groups, egcCompressedX, atom1) == 0)
+ if (getGroupType(&mtop->groups, egcCompressedX, atom2) == 0)
{
- if (getGroupType(&mtop->groups, egcCompressedX, atom2) == 0)
- {
- tng_molecule_bond_add(tng, mol, atom1,
- atom2, &tngBond);
- }
- if (getGroupType(&mtop->groups, egcCompressedX, atom3) == 0)
- {
- tng_molecule_bond_add(tng, mol, atom1,
- atom3, &tngBond);
- }
+ tng_molecule_bond_add(tng, mol, atom1,
+ atom2, &tngBond);
+ }
+ if (getGroupType(&mtop->groups, egcCompressedX, atom3) == 0)
+ {
+ tng_molecule_bond_add(tng, mol, atom1,
+ atom3, &tngBond);
}
- l += 4;
}
}
}
#include <algorithm>
#include <vector>
+#include "gromacs/compat/make_unique.h"
#include "gromacs/fileio/filetypes.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/gmxfio-xdr.h"
}
}
-static void do_ilist(t_fileio *fio, t_ilist *ilist, gmx_bool bRead)
+static void do_ilist(t_fileio *fio, InteractionList *ilist, gmx_bool bRead)
{
- gmx_fio_do_int(fio, ilist->nr);
+ int nr = ilist->size();
+ gmx_fio_do_int(fio, nr);
if (bRead)
{
- snew(ilist->iatoms, ilist->nr);
+ ilist->iatoms.resize(nr);
}
- gmx_fio_ndo_int(fio, ilist->iatoms, ilist->nr);
+ gmx_fio_ndo_int(fio, ilist->iatoms.data(), ilist->size());
}
static void do_ffparams(t_fileio *fio, gmx_ffparams_t *ffparams,
}
}
-static void add_settle_atoms(t_ilist *ilist)
+static void add_settle_atoms(InteractionList *ilist)
{
int i;
/* Settle used to only store the first atom: add the other two */
- srenew(ilist->iatoms, 2*ilist->nr);
- for (i = ilist->nr/2-1; i >= 0; i--)
+ ilist->iatoms.resize(2*ilist->size());
+ for (i = ilist->size()/4 - 1; i >= 0; i--)
{
ilist->iatoms[4*i+0] = ilist->iatoms[2*i+0];
ilist->iatoms[4*i+1] = ilist->iatoms[2*i+1];
ilist->iatoms[4*i+2] = ilist->iatoms[2*i+1] + 1;
ilist->iatoms[4*i+3] = ilist->iatoms[2*i+1] + 2;
}
- ilist->nr = 2*ilist->nr;
}
-static void do_ilists(t_fileio *fio, t_ilist *ilist, gmx_bool bRead,
+static void do_ilists(t_fileio *fio, InteractionList *ilist, gmx_bool bRead,
int file_version)
{
int j;
}
if (bClear)
{
- ilist[j].nr = 0;
- ilist[j].iatoms = nullptr;
+ ilist[j].iatoms.clear();
}
else
{
do_ilist(fio, &ilist[j], bRead);
- if (file_version < 78 && j == F_SETTLE && ilist[j].nr > 0)
+ if (file_version < 78 && j == F_SETTLE && ilist[j].size() > 0)
{
add_settle_atoms(&ilist[j]);
}
static void set_disres_npair(gmx_mtop_t *mtop)
{
- t_iparams *ip;
- gmx_mtop_ilistloop_t iloop;
- const t_ilist *ilist, *il;
- int nmol, i, npair;
- t_iatom *a;
+ t_iparams *ip;
+ gmx_mtop_ilistloop_t iloop;
+ const InteractionList *ilist;
+ int nmol;
ip = mtop->ffparams.iparams;
iloop = gmx_mtop_ilistloop_init(mtop);
while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
{
- il = &ilist[F_DISRES];
+ const InteractionList &il = ilist[F_DISRES];
- if (il->nr > 0)
+ if (il.size() > 0)
{
- a = il->iatoms;
- npair = 0;
- for (i = 0; i < il->nr; i += 3)
+ gmx::ArrayRef<const int> a = il.iatoms;
+ int npair = 0;
+ for (int i = 0; i < il.size(); i += 3)
{
npair++;
- if (i+3 == il->nr || ip[a[i]].disres.label != ip[a[i+3]].disres.label)
+ if (i+3 == il.size() || ip[a[i]].disres.label != ip[a[i+3]].disres.label)
{
ip[a[i]].disres.npair = npair;
npair = 0;
{
if (bRead)
{
- snew(mtop->intermolecular_ilist, F_NRE);
+ mtop->intermolecular_ilist = gmx::compat::make_unique<InteractionLists>();
}
- do_ilists(fio, mtop->intermolecular_ilist, bRead, file_version);
+ do_ilists(fio, mtop->intermolecular_ilist->data(), bRead, file_version);
}
}
else
{
if (fileVersion < 57)
{
- if (mtop->moltype[0].ilist[F_DISRES].nr > 0)
+ if (mtop->moltype[0].ilist[F_DISRES].size() > 0)
{
ir->eDisre = edrSimple;
}
for (i = 0; (i < asize(harm_func)); i++)
{
ft = harm_func[i];
- nh += mt->ilist[ft].nr/(interaction_function[ft].nratoms+1);
+ nh += mt->ilist[ft].size()/(interaction_function[ft].nratoms+1);
}
return nh;
}
#include <cmath>
#include <cstring>
+#include "gromacs/compat/make_unique.h"
#include "gromacs/gmxpreprocess/gpp_atomtype.h"
#include "gromacs/gmxpreprocess/topio.h"
#include "gromacs/gmxpreprocess/toputil.h"
return type;
}
-static void append_interaction(t_ilist *ilist,
+static void append_interaction(InteractionList *ilist,
int type, int nral, const int a[MAXATOMLIST])
{
- int i, where1;
-
- where1 = ilist->nr;
- ilist->nr += nral+1;
-
- ilist->iatoms[where1++] = type;
- for (i = 0; (i < nral); i++)
+ ilist->iatoms.push_back(type);
+ for (int i = 0; (i < nral); i++)
{
- ilist->iatoms[where1++] = a[i];
+ ilist->iatoms.push_back(a[i]);
}
}
static void enter_function(t_params *p, t_functype ftype, int comb, real reppow,
- gmx_ffparams_t *ffparams, t_ilist *il,
+ gmx_ffparams_t *ffparams, InteractionList *il,
int *maxtypes,
bool bNB, bool bAppend)
{
- int k, type, nr, nral, delta, start;
+ int k, type, nr, nral, start;
start = ffparams->ntypes;
nr = p->nr;
{
assert(il);
nral = NRAL(ftype);
- delta = nr*(nral+1);
- srenew(il->iatoms, il->nr+delta);
append_interaction(il, type, nral, p->param[k].a);
}
}
molt = &mtop->moltype[mt];
for (i = 0; (i < F_NRE); i++)
{
- molt->ilist[i].nr = 0;
- molt->ilist[i].iatoms = nullptr;
+ molt->ilist[i].iatoms.clear();
plist = mi[mt].plist;
if (intermolecular_interactions != nullptr)
{
/* Process the intermolecular interaction list */
- snew(mtop->intermolecular_ilist, F_NRE);
+ mtop->intermolecular_ilist = gmx::compat::make_unique<InteractionLists>();
for (i = 0; (i < F_NRE); i++)
{
- mtop->intermolecular_ilist[i].nr = 0;
- mtop->intermolecular_ilist[i].iatoms = nullptr;
+ (*mtop->intermolecular_ilist)[i].iatoms.clear();
plist = intermolecular_interactions->plist;
else
{
enter_function(&(plist[i]), static_cast<t_functype>(i), comb, reppow,
- ffp, &mtop->intermolecular_ilist[i],
+ ffp, &(*mtop->intermolecular_ilist)[i],
&maxtypes, FALSE, FALSE);
mtop->bIntermolecularInteractions = TRUE;
if (!mtop->bIntermolecularInteractions)
{
- sfree(mtop->intermolecular_ilist);
+ mtop->intermolecular_ilist.reset(nullptr);
}
}
const gmx_moltype_t *w_moltype = nullptr;
for (const gmx_moltype_t &moltype : mtop->moltype)
{
- const t_atom *atom = moltype.atoms.atom;
- const t_ilist *ilist = moltype.ilist;
- const t_ilist *ilc = &ilist[F_CONSTR];
- const t_ilist *ils = &ilist[F_SETTLE];
+ const t_atom *atom = moltype.atoms.atom;
+ const InteractionList *ilist = moltype.ilist;
+ const InteractionList *ilc = &ilist[F_CONSTR];
+ const InteractionList *ils = &ilist[F_SETTLE];
for (ftype = 0; ftype < F_NRE; ftype++)
{
if (!(ftype == F_BONDS || ftype == F_G96BONDS || ftype == F_HARMONIC))
continue;
}
- const t_ilist *ilb = &ilist[ftype];
- for (i = 0; i < ilb->nr; i += 3)
+ const InteractionList *ilb = &ilist[ftype];
+ for (i = 0; i < ilb->size(); i += 3)
{
fc = ip[ilb->iatoms[i]].harmonic.krA;
re = ip[ilb->iatoms[i]].harmonic.rA;
if (period2 < limit2)
{
bFound = false;
- for (j = 0; j < ilc->nr; j += 3)
+ for (j = 0; j < ilc->size(); j += 3)
{
if ((ilc->iatoms[j+1] == a1 && ilc->iatoms[j+2] == a2) ||
(ilc->iatoms[j+1] == a2 && ilc->iatoms[j+2] == a1))
bFound = true;
}
}
- for (j = 0; j < ils->nr; j += 4)
+ for (j = 0; j < ils->size(); j += 4)
{
if ((a1 == ils->iatoms[j+1] || a1 == ils->iatoms[j+2] || a1 == ils->iatoms[j+3]) &&
(a2 == ils->iatoms[j+1] || a2 == ils->iatoms[j+2] || a2 == ils->iatoms[j+3]))
(interaction_function[ftype].flags & IF_CONSTRAINT) ||
ftype == F_SETTLE)
{
- const t_ilist *il = &molt->ilist[ftype];
- int nral = NRAL(ftype);
+ const InteractionList *il = &molt->ilist[ftype];
+ const int nral = NRAL(ftype);
- for (int i = 0; i < il->nr; i += 1 + nral)
+ for (int i = 0; i < il->size(); i += 1 + nral)
{
for (int j = 0; j < nral; j++)
{
const t_iparams *iparams,
real massFactorThreshold)
{
- if (molt->ilist[F_CONSTR].nr == 0 && molt->ilist[F_CONSTRNC].nr == 0)
+ if (molt->ilist[F_CONSTR].size() == 0 &&
+ molt->ilist[F_CONSTRNC].size() == 0)
{
return false;
}
const t_atom * atom = molt->atoms.atom;
t_blocka atomToConstraints =
- gmx::make_at2con(molt->atoms.nr,
- molt->ilist, iparams,
+ gmx::make_at2con(*molt, iparams,
gmx::FlexibleConstraintTreatment::Exclude);
bool haveDecoupledMode = false;
{
if (interaction_function[ftype].flags & IF_ATYPE)
{
- const int nral = NRAL(ftype);
- const t_ilist *il = &molt->ilist[ftype];
- for (int i = 0; i < il->nr; i += 1 + nral)
+ const int nral = NRAL(ftype);
+ const InteractionList *il = &molt->ilist[ftype];
+ for (int i = 0; i < il->size(); i += 1 + nral)
{
/* Here we check for the mass difference between the atoms
* at both ends of the angle, that the atoms at the ends have
const gmx_groups_t *groups;
pull_params_t *pull;
int natoms, ai, aj, i, j, d, g, imin, jmin;
- t_iatom *ia;
int *nrdf2, *na_vcm, na_tot;
double *nrdf_tc, *nrdf_vcm, nrdf_uc, *nrdf_vcm_sub;
ivec *dof_vcm;
{
for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
{
- ia = molt.ilist[ftype].iatoms;
- for (i = 0; i < molt.ilist[ftype].nr; )
+ gmx::ArrayRef<const int> ia = molt.ilist[ftype].iatoms;
+ for (i = 0; i < molt.ilist[ftype].size(); )
{
/* Subtract degrees of freedom for the constraints,
* if the particles still have degrees of freedom left.
* contribute to the constraints the degrees of freedom do not
* change.
*/
- ai = as + ia[1];
- aj = as + ia[2];
- if (((atom[ia[1]].ptype == eptNucleus) ||
- (atom[ia[1]].ptype == eptAtom)) &&
- ((atom[ia[2]].ptype == eptNucleus) ||
- (atom[ia[2]].ptype == eptAtom)))
+ ai = as + ia[i + 1];
+ aj = as + ia[i + 2];
+ if (((atom[ia[i + 1]].ptype == eptNucleus) ||
+ (atom[ia[i + 1]].ptype == eptAtom)) &&
+ ((atom[ia[i + 2]].ptype == eptNucleus) ||
+ (atom[ia[i + 2]].ptype == eptAtom)))
{
if (nrdf2[ai] > 0)
{
nrdf_vcm[getGroupType(groups, egcVCM, ai)] -= 0.5*imin;
nrdf_vcm[getGroupType(groups, egcVCM, aj)] -= 0.5*jmin;
}
- ia += interaction_function[ftype].nratoms+1;
i += interaction_function[ftype].nratoms+1;
}
}
- ia = molt.ilist[F_SETTLE].iatoms;
- for (i = 0; i < molt.ilist[F_SETTLE].nr; )
+ gmx::ArrayRef<const int> ia = molt.ilist[F_SETTLE].iatoms;
+ for (i = 0; i < molt.ilist[F_SETTLE].size(); )
{
/* Subtract 1 dof from every atom in the SETTLE */
for (j = 0; j < 3; j++)
{
- ai = as + ia[1+j];
+ ai = as + ia[i + 1 + j];
imin = std::min(2, nrdf2[ai]);
nrdf2[ai] -= imin;
nrdf_tc [getGroupType(groups, egcTC, ai)] -= 0.5*imin;
nrdf_vcm[getGroupType(groups, egcVCM, ai)] -= 0.5*imin;
}
- ia += 4;
i += 4;
}
as += molt.atoms.nr;
bool posres_only,
ivec AbsRef)
{
- int d, g, i;
- gmx_mtop_ilistloop_t iloop;
- const t_ilist *ilist;
- int nmol;
- t_iparams *pr;
+ int d, g, i;
+ gmx_mtop_ilistloop_t iloop;
+ const InteractionList *ilist;
+ int nmol;
+ t_iparams *pr;
clear_ivec(AbsRef);
if (nmol > 0 &&
(AbsRef[XX] == 0 || AbsRef[YY] == 0 || AbsRef[ZZ] == 0))
{
- for (i = 0; i < ilist[F_POSRES].nr; i += 2)
+ for (i = 0; i < ilist[F_POSRES].size(); i += 2)
{
pr = &sys->ffparams.iparams[ilist[F_POSRES].iatoms[i]];
for (d = 0; d < DIM; d++)
}
}
}
- for (i = 0; i < ilist[F_FBPOSRES].nr; i += 2)
+ for (i = 0; i < ilist[F_FBPOSRES].size(); i += 2)
{
/* Check for flat-bottom posres */
pr = &sys->ffparams.iparams[ilist[F_FBPOSRES].iatoms[i]];
*/
int ftype_connbond = 0;
int ind_connbond = 0;
- if (molt->ilist[F_CONNBONDS].nr != 0)
+ if (molt->ilist[F_CONNBONDS].size() != 0)
{
fprintf(stderr, "nr. of CONNBONDS present already: %d\n",
- molt->ilist[F_CONNBONDS].nr/3);
+ molt->ilist[F_CONNBONDS].size()/3);
ftype_connbond = molt->ilist[F_CONNBONDS].iatoms[0];
- ind_connbond = molt->ilist[F_CONNBONDS].nr;
+ ind_connbond = molt->ilist[F_CONNBONDS].size();
}
/* now we delete all bonded interactions, except the ones describing
* a chemical bond. These are converted to CONNBONDS
}
int nratoms = interaction_function[ftype].nratoms;
int j = 0;
- while (j < molt->ilist[ftype].nr)
+ while (j < molt->ilist[ftype].size())
{
bool bexcl;
*/
if (bexcl && IS_CHEMBOND(ftype))
{
- srenew(molt->ilist[F_CONNBONDS].iatoms, ind_connbond + 3);
- molt->ilist[F_CONNBONDS].nr += 3;
- molt->ilist[F_CONNBONDS].iatoms[ind_connbond++] = ftype_connbond;
- molt->ilist[F_CONNBONDS].iatoms[ind_connbond++] = a1;
- molt->ilist[F_CONNBONDS].iatoms[ind_connbond++] = a2;
+ InteractionList &ilist = molt->ilist[F_CONNBONDS];
+ ilist.iatoms.resize(ind_connbond + 3);
+ ilist.iatoms[ind_connbond++] = ftype_connbond;
+ ilist.iatoms[ind_connbond++] = a1;
+ ilist.iatoms[ind_connbond++] = a2;
}
}
else
/* since the interaction involves QM atoms, these should be
* removed from the MM ilist
*/
- molt->ilist[ftype].nr -= (nratoms+1);
- for (int l = j; l < molt->ilist[ftype].nr; l++)
+ InteractionList &ilist = molt->ilist[ftype];
+ for (int k = j; k < ilist.size() - (nratoms + 1); k++)
{
- molt->ilist[ftype].iatoms[l] = molt->ilist[ftype].iatoms[l+(nratoms+1)];
+ ilist.iatoms[k] = ilist.iatoms[k + (nratoms + 1)];
}
+ ilist.iatoms.resize(ilist.size() - (nratoms + 1));
}
else
{
if (IS_CHEMBOND(i))
{
int j = 0;
- while (j < molt->ilist[i].nr)
+ while (j < molt->ilist[i].size())
{
int a1 = molt->ilist[i].iatoms[j+1];
int a2 = molt->ilist[i].iatoms[j+2];
{
int nratoms = interaction_function[i].nratoms;
int j = 0;
- while (j < molt->ilist[i].nr)
+ while (j < molt->ilist[i].size())
{
int a1 = molt->ilist[i].iatoms[j+1];
int a2 = molt->ilist[i].iatoms[j+2];
/* since the interaction involves QM atoms, these should be
* removed from the MM ilist
*/
- molt->ilist[i].nr -= (nratoms+1);
- for (int k = j; k < molt->ilist[i].nr; k++)
+ InteractionList &ilist = molt->ilist[i];
+ for (int k = j; k < ilist.size() - (nratoms + 1); k++)
{
- molt->ilist[i].iatoms[k] = molt->ilist[i].iatoms[k+(nratoms+1)];
+ ilist.iatoms[k] = ilist.iatoms[k + (nratoms + 1)];
}
+ ilist.iatoms.resize(ilist.size() - (nratoms + 1));
}
else
{
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/topology/ifunc.h"
#include "gromacs/topology/topology.h"
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/gmxassert.h"
void set_vsites_ptype(bool bVerbose, gmx_moltype_t *molt)
{
int ftype, i;
- int nra, nrd;
- t_ilist *il;
- t_iatom *ia, avsite;
if (bVerbose)
{
}
for (ftype = 0; ftype < F_NRE; ftype++)
{
- il = &molt->ilist[ftype];
+ InteractionList *il = &molt->ilist[ftype];
if (interaction_function[ftype].flags & IF_VSITE)
{
- nra = interaction_function[ftype].nratoms;
- nrd = il->nr;
- ia = il->iatoms;
+ const int nra = interaction_function[ftype].nratoms;
+ const int nrd = il->size();
+ gmx::ArrayRef<const int> ia = il->iatoms;
if (debug && nrd)
{
for (i = 0; (i < nrd); )
{
/* The virtual site */
- avsite = ia[1];
+ int avsite = ia[i + 1];
molt->atoms.atom[avsite].ptype = eptVSite;
i += nra+1;
- ia += nra+1;
}
}
}
const gmx_multisim_t *ms,
t_fcdata *fcd, t_state *state, gmx_bool bIsREMD)
{
- int fa, nmol, npair, np;
- t_disresdata *dd;
- history_t *hist;
- gmx_mtop_ilistloop_t iloop;
- const t_ilist *il;
- char *ptr;
- int type_min, type_max;
+ int fa, nmol, npair, np;
+ t_disresdata *dd;
+ history_t *hist;
+ gmx_mtop_ilistloop_t iloop;
+ const InteractionList *il;
+ char *ptr;
+ int type_min, type_max;
dd = &(fcd->disres);
iloop = gmx_mtop_ilistloop_init(mtop);
while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
{
- if (nmol > 1 && il[F_DISRES].nr > 0 && ir->eDisre != edrEnsemble)
+ if (nmol > 1 && il[F_DISRES].size() > 0 && ir->eDisre != edrEnsemble)
{
gmx_fatal(FARGS, "NMR distance restraints with multiple copies of the same molecule are currently only supported with ensemble averaging. If you just want to restrain distances between atom pairs using a flat-bottomed potential, use a restraint potential (bonds type 10) instead.");
}
np = 0;
- for (fa = 0; fa < il[F_DISRES].nr; fa += 3)
+ for (fa = 0; fa < il[F_DISRES].size(); fa += 3)
{
int type;
od->eig = nullptr;
od->v = nullptr;
- int *nr_ex = nullptr;
- int typeMin = INT_MAX;
- int typeMax = 0;
- gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(mtop);
- const t_ilist *il;
- int nmol;
+ int *nr_ex = nullptr;
+ int typeMin = INT_MAX;
+ int typeMax = 0;
+ gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(mtop);
+ const InteractionList *il;
+ int nmol;
while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
{
if (nmol > 1)
gmx_fatal(FARGS, "Found %d copies of a molecule with orientation restrains while the current code only supports a single copy. If you want to ensemble average, run multiple copies of the system using the multi-sim feature of mdrun.", nmol);
}
- for (int i = 0; i < il[F_ORIRES].nr; i += 3)
+ for (int i = 0; i < il[F_ORIRES].size(); i += 3)
{
int type = il[F_ORIRES].iatoms[i];
int ex = mtop->ffparams.iparams[type].orires.ex;
#include <cstring>
+#include "gromacs/compat/make_unique.h"
#include "gromacs/gmxlib/network.h"
#include "gromacs/math/vec.h"
#include "gromacs/mdlib/mdrun.h"
}
}
-static void bc_ilists(const t_commrec *cr, t_ilist *ilist)
+static void bc_ilists(const t_commrec *cr, InteractionList *ilist)
{
int ftype;
{
for (ftype = 0; ftype < F_NRE; ftype++)
{
- if (ilist[ftype].nr > 0)
+ if (ilist[ftype].size() > 0)
{
block_bc(cr, ftype);
- block_bc(cr, ilist[ftype].nr);
- nblock_bc(cr, ilist[ftype].nr, ilist[ftype].iatoms);
+ int nr = ilist[ftype].size();
+ block_bc(cr, nr);
+ nblock_bc(cr, nr, ilist[ftype].iatoms.data());
}
}
ftype = -1;
{
for (ftype = 0; ftype < F_NRE; ftype++)
{
- ilist[ftype].nr = 0;
+ ilist[ftype].iatoms.clear();
}
do
{
block_bc(cr, ftype);
if (ftype >= 0)
{
- block_bc(cr, ilist[ftype].nr);
- snew_bc(cr, ilist[ftype].iatoms, ilist[ftype].nr);
- nblock_bc(cr, ilist[ftype].nr, ilist[ftype].iatoms);
+ int nr;
+ block_bc(cr, nr);
+ ilist[ftype].iatoms.resize(nr);
+ nblock_bc(cr, nr, ilist[ftype].iatoms.data());
}
}
while (ftype >= 0);
block_bc(cr, mtop->bIntermolecularInteractions);
if (mtop->bIntermolecularInteractions)
{
- snew_bc(cr, mtop->intermolecular_ilist, F_NRE);
- bc_ilists(cr, mtop->intermolecular_ilist);
+ mtop->intermolecular_ilist = gmx::compat::make_unique<InteractionLists>();
+ bc_ilists(cr, mtop->intermolecular_ilist->data());
}
int nmolblock = mtop->molblock.size();
int *n_nonlin_vsite)
{
int ft, i;
- const t_ilist *il;
*n_nonlin_vsite = 0;
{
if (IS_VSITE(ft))
{
- il = &moltype->ilist[ft];
+ const InteractionList *il = &moltype->ilist[ft];
- for (i = 0; i < il->nr; i += 1+NRAL(ft))
+ for (i = 0; i < il->size(); i += 1+NRAL(ft))
{
const t_iparams *ip;
real inv_mass, coeff, m_aj;
verletbuf_atomtype_t *att;
int natt;
int ft, i, a1, a2, a3, a;
- const t_ilist *il;
const t_iparams *ip;
atom_nonbonded_kinetic_prop_t *prop;
real *vsite_m;
for (ft = F_CONSTR; ft <= F_CONSTRNC; ft++)
{
- il = &moltype.ilist[ft];
+ const InteractionList *il = &moltype.ilist[ft];
- for (i = 0; i < il->nr; i += 1+NRAL(ft))
+ for (i = 0; i < il->size(); i += 1+NRAL(ft))
{
ip = &mtop->ffparams.iparams[il->iatoms[i]];
a1 = il->iatoms[i+1];
}
}
- il = &moltype.ilist[F_SETTLE];
+ const InteractionList *il = &moltype.ilist[F_SETTLE];
- for (i = 0; i < il->nr; i += 1+NRAL(F_SETTLE))
+ for (i = 0; i < il->size(); i += 1+NRAL(F_SETTLE))
{
ip = &mtop->ffparams.iparams[il->iatoms[i]];
a1 = il->iatoms[i+1];
int nflexcon = 0;
//! A list of atoms to constraints for each moleculetype.
std::vector<t_blocka> at2con_mt;
- //! The size of at2settle = number of moltypes
- int n_at2settle_mt = 0;
- //! A list of atoms to settles.
- int **at2settle_mt = nullptr;
+ //! A list of atoms to settles for each moleculetype
+ std::vector < std::vector < int>> at2settle_mt;
//! Whether any SETTLES cross charge-group boundaries.
bool bInterCGsettles = false;
//! LINCS data.
}
}
-t_blocka make_at2con(int numAtoms,
- const t_ilist *ilists,
- const t_iparams *iparams,
- FlexibleConstraintTreatment flexibleConstraintTreatment)
+/*! \brief Returns a block struct to go from atoms to constraints
+ *
+ * The block struct will contain constraint indices with lower indices
+ * directly matching the order in F_CONSTR and higher indices matching
+ * the order in F_CONSTRNC offset by the number of constraints in F_CONSTR.
+ *
+ * \param[in] numAtoms The number of atoms to construct the list for
+ * \param[in] ilists The interaction lists, size F_NRE
+ * \param[in] iparams Interaction parameters, can be null when flexibleConstraintTreatment=Include
+ * \param[in] flexibleConstraintTreatment The flexible constraint treatment, see enum above
+ * \returns a block struct with all constraints for each atom
+ */
+template <typename T>
+static t_blocka
+makeAtomsToConstraintsList(int numAtoms,
+ const T *ilists,
+ const t_iparams *iparams,
+ FlexibleConstraintTreatment flexibleConstraintTreatment)
{
GMX_ASSERT(flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || iparams != nullptr, "With flexible constraint detection we need valid iparams");
for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
{
- const t_ilist &ilist = ilists[ftype];
- const int stride = 1 + NRAL(ftype);
- for (int i = 0; i < ilist.nr; i += stride)
+ const T &ilist = ilists[ftype];
+ const int stride = 1 + NRAL(ftype);
+ for (int i = 0; i < ilist.size(); i += stride)
{
if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include ||
!isConstraintFlexible(iparams, ilist.iatoms[i]))
int numConstraints = 0;
for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
{
- const t_ilist &ilist = ilists[ftype];
- const int stride = 1 + NRAL(ftype);
- for (int i = 0; i < ilist.nr; i += stride)
+ const T &ilist = ilists[ftype];
+ const int stride = 1 + NRAL(ftype);
+ for (int i = 0; i < ilist.size(); i += stride)
{
if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include ||
!isConstraintFlexible(iparams, ilist.iatoms[i]))
return at2con;
}
-int countFlexibleConstraints(const t_ilist *ilist,
- const t_iparams *iparams)
+t_blocka make_at2con(int numAtoms,
+ const t_ilist *ilist,
+ const t_iparams *iparams,
+ FlexibleConstraintTreatment flexibleConstraintTreatment)
+{
+ return makeAtomsToConstraintsList(numAtoms, ilist, iparams, flexibleConstraintTreatment);
+}
+
+t_blocka make_at2con(const gmx_moltype_t &moltype,
+ const t_iparams *iparams,
+ FlexibleConstraintTreatment flexibleConstraintTreatment)
+{
+ return makeAtomsToConstraintsList(moltype.atoms.nr, moltype.ilist, iparams, flexibleConstraintTreatment);
+}
+
+//! Return the number of flexible constraints in the \c ilist and \c iparams.
+template <typename T>
+static int
+countFlexibleConstraintsTemplate(const T *ilist,
+ const t_iparams *iparams)
{
int nflexcon = 0;
for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
{
const int numIatomsPerConstraint = 3;
- int ncon = ilist[ftype].nr / numIatomsPerConstraint;
- t_iatom *ia = ilist[ftype].iatoms;
- for (int con = 0; con < ncon; con++)
+ for (int i = 0; i < ilist[ftype].size(); i += numIatomsPerConstraint)
{
- if (iparams[ia[0]].constr.dA == 0 &&
- iparams[ia[0]].constr.dB == 0)
+ const int type = ilist[ftype].iatoms[i];
+ if (iparams[type].constr.dA == 0 &&
+ iparams[type].constr.dB == 0)
{
nflexcon++;
}
- ia += numIatomsPerConstraint;
}
}
return nflexcon;
}
-//! Returns the index of the settle to which each atom belongs.
-static int *make_at2settle(int natoms, const t_ilist *ilist)
+int countFlexibleConstraints(const t_ilist *ilist,
+ const t_iparams *iparams)
{
- int *at2s;
- int a, stride, s;
+ return countFlexibleConstraintsTemplate(ilist, iparams);
+}
- snew(at2s, natoms);
+//! Returns the index of the settle to which each atom belongs.
+static std::vector<int> make_at2settle(int natoms,
+ const InteractionList &ilist)
+{
/* Set all to no settle */
- for (a = 0; a < natoms; a++)
- {
- at2s[a] = -1;
- }
+ std::vector<int> at2s(natoms, -1);
- stride = 1 + NRAL(F_SETTLE);
+ const int stride = 1 + NRAL(F_SETTLE);
- for (s = 0; s < ilist->nr; s += stride)
+ for (int s = 0; s < ilist.size(); s += stride)
{
- at2s[ilist->iatoms[s+1]] = s/stride;
- at2s[ilist->iatoms[s+2]] = s/stride;
- at2s[ilist->iatoms[s+3]] = s/stride;
+ at2s[ilist.iatoms[s + 1]] = s/stride;
+ at2s[ilist.iatoms[s + 2]] = s/stride;
+ at2s[ilist.iatoms[s + 3]] = s/stride;
}
return at2s;
mapping.reserve(mtop.moltype.size());
for (const gmx_moltype_t &moltype : mtop.moltype)
{
- mapping.push_back(make_at2con(moltype.atoms.nr,
- moltype.ilist,
+ mapping.push_back(make_at2con(moltype,
mtop.ffparams.iparams,
flexibleConstraintTreatment));
}
for (const gmx_molblock_t &molblock : mtop.molblock)
{
- int count = countFlexibleConstraints(mtop.moltype[molblock.type].ilist,
+ int count =
+ countFlexibleConstraintsTemplate(mtop.moltype[molblock.type].ilist,
mtop.ffparams.iparams);
nflexcon += molblock.nmol*count;
}
settled = settle_init(mtop);
/* Make an atom to settle index for use in domain decomposition */
- n_at2settle_mt = mtop.moltype.size();
- snew(at2settle_mt, n_at2settle_mt);
for (size_t mt = 0; mt < mtop.moltype.size(); mt++)
{
- at2settle_mt[mt] =
- make_at2settle(mtop.moltype[mt].atoms.nr,
- &mtop.moltype[mt].ilist[F_SETTLE]);
+ at2settle_mt.emplace_back(make_at2settle(mtop.moltype[mt].atoms.nr,
+ mtop.moltype[mt].ilist[F_SETTLE]));
}
/* Allocate thread-local work arrays */
return impl_->at2con_mt;
}
-int *const* Constraints::atom2settle_moltype() const
+ArrayRef < const std::vector < int>> Constraints::atom2settle_moltype() const
{
return impl_->at2settle_mt;
}
{
const gmx_moltype_t *molt;
const t_block *cgs;
- const t_ilist *il;
int *at2cg, cg, a, ftype, i;
bool bInterCG;
{
molt = &mtop.moltype[mtop.molblock[mb].type];
- if (molt->ilist[F_CONSTR].nr > 0 ||
- molt->ilist[F_CONSTRNC].nr > 0 ||
- molt->ilist[F_SETTLE].nr > 0)
+ if (molt->ilist[F_CONSTR].size() > 0 ||
+ molt->ilist[F_CONSTRNC].size() > 0 ||
+ molt->ilist[F_SETTLE].size() > 0)
{
cgs = &molt->cgs;
snew(at2cg, molt->atoms.nr);
for (ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
{
- il = &molt->ilist[ftype];
- for (i = 0; i < il->nr && !bInterCG; i += 1+NRAL(ftype))
+ const InteractionList &il = molt->ilist[ftype];
+ for (i = 0; i < il.size() && !bInterCG; i += 1+NRAL(ftype))
{
- if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]])
+ if (at2cg[il.iatoms[i + 1]] != at2cg[il.iatoms[i + 2]])
{
bInterCG = TRUE;
}
{
const gmx_moltype_t *molt;
const t_block *cgs;
- const t_ilist *il;
int *at2cg, cg, a, ftype, i;
bool bInterCG;
{
molt = &mtop.moltype[mtop.molblock[mb].type];
- if (molt->ilist[F_SETTLE].nr > 0)
+ if (molt->ilist[F_SETTLE].size() > 0)
{
cgs = &molt->cgs;
snew(at2cg, molt->atoms.nr);
for (ftype = F_SETTLE; ftype <= F_SETTLE; ftype++)
{
- il = &molt->ilist[ftype];
- for (i = 0; i < il->nr && !bInterCG; i += 1+NRAL(F_SETTLE))
+ const InteractionList &il = molt->ilist[ftype];
+ for (i = 0; i < il.size() && !bInterCG; i += 1+NRAL(F_SETTLE))
{
- if (at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+2]] ||
- at2cg[il->iatoms[i+1]] != at2cg[il->iatoms[i+3]])
+ if (at2cg[il.iatoms[i + 1]] != at2cg[il.iatoms[i + 2]] ||
+ at2cg[il.iatoms[i + 1]] != at2cg[il.iatoms[i + 3]])
{
bInterCG = TRUE;
}
struct gmx_edsam;
struct gmx_localtop_t;
+struct gmx_moltype_t;
struct gmx_mtop_t;
struct gmx_multisim_t;
struct gmx_wallcycle;
//! Getter for use by domain decomposition.
const ArrayRef<const t_blocka> atom2constraints_moltype() const;
//! Getter for use by domain decomposition.
- int *const* atom2settle_moltype() const;
+ ArrayRef < const std::vector < int>> atom2settle_moltype() const;
/*! \brief Return the data for reduction for determining
* constraint RMS relative deviations, or an empty ArrayRef
FlexibleConstraintTreatment
flexibleConstraintTreatment(bool haveDynamicsIntegrator);
+/*! \brief Returns a block struct to go from atoms to constraints
+ *
+ * The block struct will contain constraint indices with lower indices
+ * directly matching the order in F_CONSTR and higher indices matching
+ * the order in F_CONSTRNC offset by the number of constraints in F_CONSTR.
+ *
+ * \param[in] moltype The molecule data
+ * \param[in] iparams Interaction parameters, can be null when flexibleConstraintTreatment=Include
+ * \param[in] flexibleConstraintTreatment The flexible constraint treatment, see enum above
+ * \returns a block struct with all constraints for each atom
+ */
+t_blocka make_at2con(const gmx_moltype_t &moltype,
+ const t_iparams *iparams,
+ FlexibleConstraintTreatment flexibleConstraintTreatment);
+
/*! \brief Returns a block struct to go from atoms to constraints
*
* The block struct will contain constraint indices with lower indices
int countFlexibleConstraints(const t_ilist *ilist,
const t_iparams *iparams);
-/*! \brief Macro for getting the constraint iatoms for a constraint number con
+/*! \brief Returns the constraint iatoms for a constraint number con
* which comes from a list where F_CONSTR and F_CONSTRNC constraints
* are concatenated. */
-#define constr_iatomptr(nconstr, iatom_constr, iatom_constrnc, con) ((con) < (nconstr) ? (iatom_constr)+(con)*3 : (iatom_constrnc)+((con)-(nconstr))*3)
+inline const int *
+constr_iatomptr(gmx::ArrayRef<const int> iatom_constr,
+ gmx::ArrayRef<const int> iatom_constrnc,
+ int con)
+{
+ if (con*3 < iatom_constr.size())
+ {
+ return iatom_constr.data() + con*3;
+ }
+ else
+ {
+ return iatom_constrnc.data() + con*3 - iatom_constr.size();
+ }
+};
/*! \brief Returns whether there are inter charge group constraints */
bool inter_charge_group_constraints(const gmx_mtop_t &mtop);
//! Recursing function to help find all adjacent constraints.
static void constr_recur(const t_blocka *at2con,
- const t_ilist *ilist, const t_iparams *iparams,
+ const InteractionList *ilist, const t_iparams *iparams,
gmx_bool bTopB,
int at, int depth, int nc, int *path,
real r0, real r1, real *r2max,
int *count)
{
- int ncon1;
- t_iatom *ia1, *ia2;
int c, con, a1;
gmx_bool bUse;
- t_iatom *ia;
real len, rn0, rn1;
(*count)++;
- ncon1 = ilist[F_CONSTR].nr/3;
- ia1 = ilist[F_CONSTR].iatoms;
- ia2 = ilist[F_CONSTRNC].iatoms;
+ gmx::ArrayRef<const int> ia1 = ilist[F_CONSTR].iatoms;
+ gmx::ArrayRef<const int> ia2 = ilist[F_CONSTRNC].iatoms;
/* Loop over all constraints connected to this atom */
for (c = at2con->index[at]; c < at2con->index[at+1]; c++)
}
if (bUse)
{
- ia = constr_iatomptr(ncon1, ia1, ia2, con);
+ const int *ia = constr_iatomptr(ia1, ia2, con);
/* Flexible constraints currently have length 0, which is incorrect */
if (!bTopB)
{
{
fprintf(debug, " %d %5.3f",
path[a1],
- iparams[constr_iatomptr(ncon1, ia1, ia2, con)[0]].constr.dA);
+ iparams[constr_iatomptr(ia1, ia2, con)[0]].constr.dA);
}
fprintf(debug, " %d %5.3f\n", con, len);
}
t_blocka at2con;
real r0, r1, r2maxA, r2maxB, rmax, lam0, lam1;
- if (molt->ilist[F_CONSTR].nr == 0 &&
- molt->ilist[F_CONSTRNC].nr == 0)
+ if (molt->ilist[F_CONSTR].size() == 0 &&
+ molt->ilist[F_CONSTRNC].size() == 0)
{
return 0;
}
natoms = molt->atoms.nr;
- at2con = make_at2con(natoms, molt->ilist, iparams,
+ at2con = make_at2con(*molt, iparams,
flexibleConstraintTreatment(EI_DYNAMICS(ir->eI)));
snew(path, 1+ir->nProjOrder);
for (at = 0; at < 1+ir->nProjOrder; at++)
int nral;
nral = NRAL(ftype);
- for (ia = 0; ia < molt->ilist[ftype].nr; ia += 1+nral)
+ for (ia = 0; ia < molt->ilist[ftype].size(); ia += 1+nral)
{
int a;
static void count_tables(int ftype1, int ftype2, const gmx_mtop_t *mtop,
int *ncount, int **count)
{
- const t_ilist *il;
- int ftype, stride, i, j, tabnr;
+ int ftype, i, j, tabnr;
// Loop over all moleculetypes
for (const gmx_moltype_t &molt : mtop->moltype)
// If the current interaction type is one of the types whose tables we're trying to count...
if (ftype == ftype1 || ftype == ftype2)
{
- il = &molt.ilist[ftype];
- stride = 1 + NRAL(ftype);
+ const InteractionList &il = molt.ilist[ftype];
+ const int stride = 1 + NRAL(ftype);
// ... and there are actually some interactions for this type
- for (i = 0; i < il->nr; i += stride)
+ for (i = 0; i < il.size(); i += stride)
{
// Find out which table index the user wanted
- tabnr = mtop->ffparams.iparams[il->iatoms[i]].tab.table;
+ tabnr = mtop->ffparams.iparams[il.iatoms[i]].tab.table;
if (tabnr < 0)
{
gmx_fatal(FARGS, "A bonded table number is smaller than 0: %d\n", tabnr);
}
//! Finds all triangles of atoms that share constraints to a central atom.
-static int count_triangle_constraints(const t_ilist *ilist,
- const t_blocka &at2con)
+static int count_triangle_constraints(const InteractionList *ilist,
+ const t_blocka &at2con)
{
int ncon1, ncon_tot;
- int c0, a00, a01, n1, c1, a10, a11, ac1, n2, c2, a20, a21;
+ int c0, n1, c1, ac1, n2, c2;
int ncon_triangle;
- bool bTriangle;
- t_iatom *ia1, *ia2, *iap;
- ncon1 = ilist[F_CONSTR].nr/3;
- ncon_tot = ncon1 + ilist[F_CONSTRNC].nr/3;
+ ncon1 = ilist[F_CONSTR].size()/3;
+ ncon_tot = ncon1 + ilist[F_CONSTRNC].size()/3;
- ia1 = ilist[F_CONSTR].iatoms;
- ia2 = ilist[F_CONSTRNC].iatoms;
+ gmx::ArrayRef<const int> ia1 = ilist[F_CONSTR].iatoms;
+ gmx::ArrayRef<const int> ia2 = ilist[F_CONSTRNC].iatoms;
ncon_triangle = 0;
for (c0 = 0; c0 < ncon_tot; c0++)
{
- bTriangle = FALSE;
- iap = constr_iatomptr(ncon1, ia1, ia2, c0);
- a00 = iap[1];
- a01 = iap[2];
+ bool bTriangle = FALSE;
+ const int *iap = constr_iatomptr(ia1, ia2, c0);
+ const int a00 = iap[1];
+ const int a01 = iap[2];
for (n1 = at2con.index[a01]; n1 < at2con.index[a01+1]; n1++)
{
c1 = at2con.a[n1];
if (c1 != c0)
{
- iap = constr_iatomptr(ncon1, ia1, ia2, c1);
- a10 = iap[1];
- a11 = iap[2];
+ const int *iap = constr_iatomptr(ia1, ia2, c1);
+ const int a10 = iap[1];
+ const int a11 = iap[2];
if (a10 == a01)
{
ac1 = a11;
c2 = at2con.a[n2];
if (c2 != c0 && c2 != c1)
{
- iap = constr_iatomptr(ncon1, ia1, ia2, c2);
- a20 = iap[1];
- a21 = iap[2];
+ const int *iap = constr_iatomptr(ia1, ia2, c2);
+ const int a20 = iap[1];
+ const int a21 = iap[2];
if (a20 == a00 || a21 == a00)
{
bTriangle = TRUE;
}
//! Finds sequences of sequential constraints.
-static bool more_than_two_sequential_constraints(const t_ilist *ilist,
- const t_blocka &at2con)
+static bool more_than_two_sequential_constraints(const InteractionList *ilist,
+ const t_blocka &at2con)
{
- t_iatom *ia1, *ia2, *iap;
int ncon1, ncon_tot, c;
- int a1, a2;
bool bMoreThanTwoSequentialConstraints;
- ncon1 = ilist[F_CONSTR].nr/3;
- ncon_tot = ncon1 + ilist[F_CONSTRNC].nr/3;
+ ncon1 = ilist[F_CONSTR].size()/3;
+ ncon_tot = ncon1 + ilist[F_CONSTRNC].size()/3;
- ia1 = ilist[F_CONSTR].iatoms;
- ia2 = ilist[F_CONSTRNC].iatoms;
+ gmx::ArrayRef<const int> ia1 = ilist[F_CONSTR].iatoms;
+ gmx::ArrayRef<const int> ia2 = ilist[F_CONSTRNC].iatoms;
bMoreThanTwoSequentialConstraints = FALSE;
for (c = 0; c < ncon_tot && !bMoreThanTwoSequentialConstraints; c++)
{
- iap = constr_iatomptr(ncon1, ia1, ia2, c);
- a1 = iap[1];
- a2 = iap[2];
+ const int *iap = constr_iatomptr(ia1, ia2, c);
+ const int a1 = iap[1];
+ const int a2 = iap[2];
/* Check if this constraint has constraints connected at both atoms */
if (at2con.index[a1+1] - at2con.index[a1] > 1 &&
at2con.index[a2+1] - at2con.index[a2] > 1)
{
for (j = 0; j < F_LJ; j++)
{
- mtop->moltype[i].ilist[j].nr = 0;
+ mtop->moltype[i].ilist[j].iatoms.clear();
}
for (j = F_POSRES; j <= F_VSITEN; j++)
{
- mtop->moltype[i].ilist[j].nr = 0;
+ mtop->moltype[i].ilist[j].iatoms.clear();
}
}
}
nd_c = NRAL(ftype) - 1;
break;
}
- nbonds = molb.nmol*molt->ilist[ftype].nr/(1 + NRAL(ftype));
+ nbonds = molb.nmol*molt->ilist[ftype].size()/(1 + NRAL(ftype));
ndtot_c += nbonds*nd_c;
ndtot_simd += nbonds*nd_simd;
}
settledata *settle_init(const gmx_mtop_t &mtop)
{
/* Check that we have only one settle type */
- int settle_type = -1;
- gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(mtop);
- const t_ilist *ilist;
- int nmol;
- const int nral1 = 1 + NRAL(F_SETTLE);
+ int settle_type = -1;
+ gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(mtop);
+ const InteractionList *ilist;
+ int nmol;
+ const int nral1 = 1 + NRAL(F_SETTLE);
while (gmx_mtop_ilistloop_next(iloop, &ilist, &nmol))
{
- for (int i = 0; i < ilist[F_SETTLE].nr; i += nral1)
+ for (int i = 0; i < ilist[F_SETTLE].size(); i += nral1)
{
if (settle_type == -1)
{
int aS, aN = 0; /* Shell and nucleus */
int bondtypes[] = { F_BONDS, F_HARMONIC, F_CUBICBONDS, F_POLARIZATION, F_ANHARM_POL, F_WATER_POL };
#define NBT asize(bondtypes)
- t_iatom *ia;
gmx_mtop_atomloop_all_t aloop;
const gmx_ffparams_t *ffparams;
{
for (j = 0; (j < NBT); j++)
{
- ia = molt->ilist[bondtypes[j]].iatoms;
- for (i = 0; (i < molt->ilist[bondtypes[j]].nr); )
+ const int *ia = molt->ilist[bondtypes[j]].iatoms.data();
+ for (i = 0; (i < molt->ilist[bondtypes[j]].size()); )
{
type = ia[0];
ftype = ffparams->functype[type];
}
else
{
- /* Pass NULL iso fplog to avoid graph prints for each molecule type */
- mk_graph_ilist(nullptr, moltype.ilist,
- 0, moltype.atoms.nr, FALSE, FALSE, graph);
+ mk_graph_moltype(moltype, graph);
for (mol = 0; mol < molb.nmol; mol++)
{
mtop.moltype.resize(1);
mtop.molblock.resize(1);
mtop.molblock[0].type = 0;
- const int settleStride = 1 + atomsPerSettle;
- int *iatoms;
- snew(iatoms, numSettles*settleStride);
+ std::vector<int> &iatoms = mtop.moltype[0].ilist[F_SETTLE].iatoms;
for (int i = 0; i < numSettles; ++i)
{
- iatoms[i*settleStride + 0] = settleType;
- iatoms[i*settleStride + 1] = i*atomsPerSettle + 0;
- iatoms[i*settleStride + 2] = i*atomsPerSettle + 1;
- iatoms[i*settleStride + 3] = i*atomsPerSettle + 2;
+ iatoms.push_back(settleType);
+ iatoms.push_back(i*atomsPerSettle + 0);
+ iatoms.push_back(i*atomsPerSettle + 1);
+ iatoms.push_back(i*atomsPerSettle + 2);
}
- mtop.moltype[0].ilist[F_SETTLE].iatoms = iatoms;
- mtop.moltype[0].ilist[F_SETTLE].nr = numSettles*settleStride;
// Set up the SETTLE parameters.
mtop.ffparams.ntypes = 1;
mdatoms.homenr = numSettles * atomsPerSettle;
// Finally make the settle data structures
- settledata *settled = settle_init(mtop);
- settle_set_constraints(settled, &mtop.moltype[0].ilist[F_SETTLE], mdatoms);
+ settledata *settled = settle_init(mtop);
+ const t_ilist ilist = { mtop.moltype[0].ilist[F_SETTLE].size(), 0, mtop.moltype[0].ilist[F_SETTLE].iatoms.data(), 0 };
+ settle_set_constraints(settled, &ilist, mdatoms);
// Copy the original positions from the array of doubles to a vector of reals
std::vector<real> startingPositions(std::begin(g_positions), std::end(g_positions));
*
* \param[in] ilist The interaction list
*/
-static int vsiteIlistNrCount(const t_ilist *ilist)
+template <typename T>
+static int vsiteIlistNrCount(const T *ilist)
{
int nr = 0;
for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
{
- nr += ilist[ftype].nr;
+ nr += ilist[ftype].size();
}
return nr;
int atomOffset = mtop.moleculeBlockIndices[mb].globalAtomStart;
for (int mol = 0; mol < molb.nmol; mol++)
{
+ t_ilist ilist[F_NRE];
+ for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
+ {
+ ilist[ftype].nr = molt.ilist[ftype].size();
+ ilist[ftype].iatoms = const_cast<t_iatom *>(molt.ilist[ftype].iatoms.data());
+ }
+
construct_vsites(nullptr, as_rvec_array(x.data()) + atomOffset,
0.0, nullptr,
- mtop.ffparams.iparams, molt.ilist,
+ mtop.ffparams.iparams, ilist,
epbcNONE, TRUE, nullptr, nullptr);
atomOffset += molt.atoms.nr;
}
std::vector<int> a2cg = atom2cg(molt.cgs);
for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
{
- int nral = NRAL(ftype);
- const t_ilist &il = molt.ilist[ftype];
- const t_iatom *ia = il.iatoms;
- for (int i = 0; i < il.nr; i += 1 + nral)
+ const int nral = NRAL(ftype);
+ const InteractionList &il = molt.ilist[ftype];
+ for (int i = 0; i < il.size(); i += 1 + nral)
{
- int cg = a2cg[ia[1+i]];
+ int cg = a2cg[il.iatoms[1 + i]];
for (int a = 1; a < nral; a++)
{
- if (a2cg[ia[1+a]] != cg)
+ if (a2cg[il.iatoms[1 + a]] != cg)
{
n_intercg_vsite += molb.nmol;
break;
return n_intercg_vsite;
}
+template <typename T>
static VsitePbc
-get_vsite_pbc(const t_iparams *iparams, const t_ilist *ilist,
+get_vsite_pbc(const t_iparams *iparams, const T *ilist,
const t_atom *atom, const t_mdatoms *md,
const t_block &cgs)
{
{
{ // TODO remove me
int nral = NRAL(ftype);
- const t_ilist *il = &ilist[ftype];
- const t_iatom *ia = il->iatoms;
+ const T &il = ilist[ftype];
std::vector<int> &vsite_pbc_f = vsite_pbc[ftype - F_VSITE2];
- vsite_pbc_f.resize(il->nr/(1 + nral));
+ vsite_pbc_f.resize(il.size()/(1 + nral));
int i = 0;
- while (i < il->nr)
+ while (i < il.size())
{
int vsi = i/(1 + nral);
- t_iatom vsite = ia[i+1];
+ t_iatom vsite = il.iatoms[i + 1];
int cg_v = a2cg[vsite];
/* A value of -2 signals that this vsite and its contructing
* atoms are all within the same cg, so no pbc is required.
int nc3 = 0;
if (ftype == F_VSITEN)
{
- nc3 = 3*iparams[ia[i]].vsiten.n;
+ nc3 = 3*iparams[il.iatoms[i]].vsiten.n;
for (int j = 0; j < nc3; j += 3)
{
- if (a2cg[ia[i+j+2]] != cg_v)
+ if (a2cg[il.iatoms[i + j + 2]] != cg_v)
{
vsite_pbc_f[vsi] = -1;
}
{
for (int a = 1; a < nral; a++)
{
- if (a2cg[ia[i+1+a]] != cg_v)
+ if (a2cg[il.iatoms[i + 1 + a]] != cg_v)
{
vsite_pbc_f[vsi] = -1;
}
*/
vsite_pbc_f[vsi] = vsite;
}
- else if (cg_v != a2cg[ia[1+i+1]])
+ else if (cg_v != a2cg[il.iatoms[1 + i + 1]])
{
/* This vsite has a different charge group index
* than it's first constructing atom
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/topology/idef.h"
#include "gromacs/topology/ifunc.h"
+#include "gromacs/topology/topology.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/strconvert.h"
* When a non-null part array is supplied with part indices for each atom,
* edges are only added when atoms have a different part index.
*/
-static bool mk_igraph(t_graph *g, int ftype, const t_ilist *il,
+template <typename T>
+static bool mk_igraph(t_graph *g, int ftype, const T &il,
int at_start, int at_end,
const int *part)
{
- t_iatom *ia;
int i, j, np;
int end;
bool addedEdge = false;
- end = il->nr;
- ia = il->iatoms;
+ end = il.size();
i = 0;
while (i < end)
{
np = interaction_function[ftype].nratoms;
- if (np > 1 && ia[1] >= at_start && ia[1] < at_end)
+ if (np > 1 && il.iatoms[i + 1] >= at_start && il.iatoms[i + 1] < at_end)
{
- if (ia[np] >= at_end)
+ if (il.iatoms[i + np] >= at_end)
{
gmx_fatal(FARGS,
"Molecule in topology has atom numbers below and "
if (ftype == F_SETTLE)
{
/* Bond all the atoms in the settle */
- add_gbond(g, ia[1], ia[2]);
- add_gbond(g, ia[1], ia[3]);
+ add_gbond(g, il.iatoms[i + 1], il.iatoms[i + 2]);
+ add_gbond(g, il.iatoms[i + 1], il.iatoms[i + 3]);
addedEdge = true;
}
else if (part == nullptr)
/* Simply add this bond */
for (j = 1; j < np; j++)
{
- add_gbond(g, ia[j], ia[j+1]);
+ add_gbond(g, il.iatoms[i + j], il.iatoms[i + j+1]);
}
addedEdge = true;
}
/* Add this bond when it connects two unlinked parts of the graph */
for (j = 1; j < np; j++)
{
- if (part[ia[j]] != part[ia[j+1]])
+ if (part[il.iatoms[i + j]] != part[il.iatoms[i + j+1]])
{
- add_gbond(g, ia[j], ia[j+1]);
+ add_gbond(g, il.iatoms[i + j], il.iatoms[i + j+1]);
addedEdge = true;
}
}
}
}
- ia += np+1;
+
i += np+1;
}
fflush(log);
}
-static void calc_1se(t_graph *g, int ftype, const t_ilist *il,
+template <typename T>
+static void calc_1se(t_graph *g, int ftype, const T &il,
int nbond[], int at_start, int at_end)
{
int k, nratoms, end, j;
- t_iatom *ia, iaa;
- end = il->nr;
+ end = il.size();
- ia = il->iatoms;
- for (j = 0; (j < end); j += nratoms+1, ia += nratoms+1)
+ for (j = 0; (j < end); j += nratoms + 1)
{
nratoms = interaction_function[ftype].nratoms;
if (ftype == F_SETTLE)
{
- iaa = ia[1];
+ const int iaa = il.iatoms[j + 1];
if (iaa >= at_start && iaa < at_end)
{
- nbond[iaa] += 2;
- nbond[ia[2]] += 1;
- nbond[ia[3]] += 1;
- g->at_start = std::min(g->at_start, iaa);
- g->at_end = std::max(g->at_end, iaa+2+1);
+ nbond[iaa] += 2;
+ nbond[il.iatoms[j + 2]] += 1;
+ nbond[il.iatoms[j + 3]] += 1;
+ g->at_start = std::min(g->at_start, iaa);
+ g->at_end = std::max(g->at_end, iaa+2+1);
}
}
else
{
for (k = 1; (k <= nratoms); k++)
{
- iaa = ia[k];
+ const int iaa = il.iatoms[j + k];
if (iaa >= at_start && iaa < at_end)
{
g->at_start = std::min(g->at_start, iaa);
}
}
-static int calc_start_end(FILE *fplog, t_graph *g, const t_ilist il[],
+template <typename T>
+static int calc_start_end(FILE *fplog, t_graph *g, const T il[],
int at_start, int at_end,
int nbond[])
{
{
if (interaction_function[i].flags & IF_CHEMBOND)
{
- calc_1se(g, i, &il[i], nbond, at_start, at_end);
+ calc_1se(g, i, il[i], nbond, at_start, at_end);
}
}
/* Then add all the other interactions in fixed lists, but first
{
if (!(interaction_function[i].flags & IF_CHEMBOND))
{
- calc_1se(g, i, &il[i], nbond, at_start, at_end);
+ calc_1se(g, i, il[i], nbond, at_start, at_end);
}
}
return bMultiPart;
}
-void mk_graph_ilist(FILE *fplog,
- const t_ilist *ilist, int at_start, int at_end,
- gmx_bool bShakeOnly, gmx_bool bSettle,
- t_graph *g)
+template <typename T>
+static void
+mk_graph_ilist(FILE *fplog,
+ const T *ilist, int at_start, int at_end,
+ gmx_bool bShakeOnly, gmx_bool bSettle,
+ t_graph *g)
{
int *nbond;
int i, nbtot;
{
if (interaction_function[i].flags & IF_CHEMBOND)
{
- mk_igraph(g, i, &(ilist[i]), at_start, at_end, nullptr);
+ mk_igraph(g, i, ilist[i], at_start, at_end, nullptr);
}
}
if (!(interaction_function[i].flags & IF_CHEMBOND))
{
bool addedEdgeForType =
- mk_igraph(g, i, &(ilist[i]), at_start, at_end, nbond);
+ mk_igraph(g, i, ilist[i], at_start, at_end, nbond);
addedEdge = (addedEdge || addedEdgeForType);
}
}
else
{
/* This is a special thing used in splitter.c to generate shake-blocks */
- mk_igraph(g, F_CONSTR, &(ilist[F_CONSTR]), at_start, at_end, nullptr);
+ mk_igraph(g, F_CONSTR, ilist[F_CONSTR], at_start, at_end, nullptr);
if (bSettle)
{
- mk_igraph(g, F_SETTLE, &(ilist[F_SETTLE]), at_start, at_end, nullptr);
+ mk_igraph(g, F_SETTLE, ilist[F_SETTLE], at_start, at_end, nullptr);
}
}
g->nbound = 0;
}
}
+void mk_graph_moltype(const gmx_moltype_t &moltype,
+ t_graph *g)
+{
+ mk_graph_ilist(nullptr, moltype.ilist, 0, moltype.atoms.nr,
+ FALSE, FALSE,
+ g);
+}
+
t_graph *mk_graph(FILE *fplog,
const t_idef *idef, int at_start, int at_end,
gmx_bool bShakeOnly, gmx_bool bSettle)
#include "gromacs/math/vectypes.h"
#include "gromacs/utility/basedefinitions.h"
+struct InteractionList;
+struct gmx_moltype_t;
struct t_idef;
-struct t_ilist;
typedef enum {
egcolWhite, egcolGrey, egcolBlack, egcolNR
* If bSettle && bShakeOnly the settles are used too.
*/
-void mk_graph_ilist(FILE *fplog,
- const struct t_ilist *ilist, int at_start, int at_end,
- gmx_bool bShakeOnly, gmx_bool bSettle,
- t_graph *g);
-/* As mk_graph, but takes t_ilist iso t_idef and does not allocate g */
+void mk_graph_moltype(const gmx_moltype_t &moltype,
+ t_graph *g);
+/* As mk_graph, but takes gmx_moltype_t iso t_idef and does not allocate g */
void done_graph(t_graph *g);
mtop->moltype[0].atoms = top.atoms;
for (i = 0; i < F_NRE; i++)
{
- mtop->moltype[0].ilist[i] = top.idef.il[i];
+ InteractionList &ilist = mtop->moltype[0].ilist[i];
+ ilist.iatoms.resize(top.idef.il[i].nr);
+ for (int j = 0; j < top.idef.il[i].nr; j++)
+ {
+ ilist.iatoms[j] = top.idef.il[i].iatoms[j];
+ }
}
mtop->moltype[0].atoms = top.atoms;
mtop->moltype[0].cgs = top.cgs;
}
}
-void pr_ilist(FILE *fp, int indent, const char *title,
- const t_functype *functype, const t_ilist *ilist,
- gmx_bool bShowNumbers,
- gmx_bool bShowParameters, const t_iparams *iparams)
+template <typename T>
+static void
+printIlist(FILE *fp, int indent, const char *title,
+ const t_functype *functype, const T *ilist,
+ gmx_bool bShowNumbers,
+ gmx_bool bShowParameters, const t_iparams *iparams)
{
int i, j, k, type, ftype;
- t_iatom *iatoms;
- if (available(fp, ilist, indent, title) && ilist->nr > 0)
+ if (available(fp, ilist, indent, title) && ilist->size() > 0)
{
indent = pr_title(fp, indent, title);
pr_indent(fp, indent);
- fprintf(fp, "nr: %d\n", ilist->nr);
- if (ilist->nr > 0)
+ fprintf(fp, "nr: %d\n", ilist->size());
+ if (ilist->size() > 0)
{
pr_indent(fp, indent);
fprintf(fp, "iatoms:\n");
- iatoms = ilist->iatoms;
- for (i = j = 0; i < ilist->nr; )
+ for (i = j = 0; i < ilist->size(); )
{
pr_indent(fp, indent+INDENT);
- type = *(iatoms++);
+ type = ilist->iatoms[i];
ftype = functype[type];
if (bShowNumbers)
{
printf("(%s)", interaction_function[ftype].name);
for (k = 0; k < interaction_function[ftype].nratoms; k++)
{
- fprintf(fp, " %3d", *(iatoms++));
+ fprintf(fp, " %3d", ilist->iatoms[i + 1 + k]);
}
if (bShowParameters)
{
}
}
+void pr_ilist(FILE *fp, int indent, const char *title,
+ const t_functype *functype, const InteractionList *ilist,
+ gmx_bool bShowNumbers,
+ gmx_bool bShowParameters, const t_iparams *iparams)
+{
+ printIlist(fp, indent, title, functype, ilist,
+ bShowNumbers, bShowParameters, iparams);
+}
+
static void pr_cmap(FILE *fp, int indent, const char *title,
const gmx_cmap_t *cmap_grid, gmx_bool bShowNumbers)
{
for (j = 0; (j < F_NRE); j++)
{
- pr_ilist(fp, indent, interaction_function[j].longname,
- idef->functype, &idef->il[j], bShowNumbers,
- bShowParameters, idef->iparams);
+ printIlist(fp, indent, interaction_function[j].longname,
+ idef->functype, &idef->il[j], bShowNumbers,
+ bShowParameters, idef->iparams);
}
}
}
#include <stdio.h>
+#include <array>
+#include <vector>
+
#include "gromacs/math/vectypes.h"
#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/real.h"
typedef int t_functype;
-/*
+/* List of listed interactions, see description further down.
+ *
+ * TODO: Consider storing the function type as well.
+ * TODO: Consider providing per interaction access.
+ */
+struct InteractionList
+{
+ /* Returns the total number of elements in iatoms */
+ int size() const
+ {
+ return iatoms.size();
+ }
+
+ /* List of interactions, see explanation further down */
+ std::vector<int> iatoms;
+};
+
+/* List of interaction lists, one list for each interaction type
+ *
+ * TODO: Consider only including entries in use instead of all F_NRE
+ */
+typedef std::array<InteractionList, F_NRE> InteractionLists;
+
+/* Deprecated list of listed interactions.
+ *
* The nonperturbed/perturbed interactions are now separated (sorted) in the
* ilist, such that the first 0..(nr_nonperturbed-1) ones are exactly that, and
* the remaining ones from nr_nonperturbed..(nr-1) are perturbed bonded
* interactions.
*/
-typedef struct t_ilist
+struct t_ilist
{
+ /* Returns the total number of elements in iatoms */
+ int size() const
+ {
+ return nr;
+ }
+
int nr;
int nr_nonperturbed;
t_iatom *iatoms;
int nalloc;
-} t_ilist;
+};
+
+/* TODO: Replace t_ilist in gmx_localtop_t by InteractionList.
+ * The nr_nonperturbed functionality needs to be ported.
+ * Remove t_topology.
+ * Remove t_ilist and remove templating on list type
+ * in mshift.cpp, constr.cpp, vsite.cpp and domdec_topology.cpp.
+ */
/*
- * The struct t_ilist defines a list of atoms with their interactions.
+ * The structs InteractionList and t_ilist defines a list of atoms with their interactions.
* General field description:
* int nr
* the size (nr elements) of the interactions array (iatoms[]).
void pr_iparams(FILE *fp, t_functype ftype, const t_iparams *iparams);
void pr_ilist(FILE *fp, int indent, const char *title,
- const t_functype *functype, const t_ilist *ilist,
+ const t_functype *functype, const InteractionList *ilist,
gmx_bool bShowNumbers,
gmx_bool bShowParameters, const t_iparams *iparams);
void pr_ffparams(FILE *fp, int indent, const char *title,
sfree(iloop);
}
-gmx_bool gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop,
- const t_ilist **ilist_mol, int *nmol)
+gmx_bool gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop,
+ const InteractionList **ilist_mol,
+ int *nmol)
{
if (iloop == nullptr)
{
if (iloop->mblock == static_cast<int>(iloop->mtop->molblock.size()) &&
iloop->mtop->bIntermolecularInteractions)
{
- *ilist_mol = iloop->mtop->intermolecular_ilist;
+ *ilist_mol = iloop->mtop->intermolecular_ilist->data();
*nmol = 1;
return TRUE;
}
sfree(iloop);
}
-gmx_bool gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop,
- const t_ilist **ilist_mol, int *atnr_offset)
+gmx_bool gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop,
+ const InteractionList **ilist_mol,
+ int *atnr_offset)
{
if (iloop == nullptr)
if (iloop->mblock == iloop->mtop->molblock.size() &&
iloop->mtop->bIntermolecularInteractions)
{
- *ilist_mol = iloop->mtop->intermolecular_ilist;
+ *ilist_mol = iloop->mtop->intermolecular_ilist->data();
*atnr_offset = 0;
return TRUE;
}
int gmx_mtop_ftype_count(const gmx_mtop_t *mtop, int ftype)
{
- gmx_mtop_ilistloop_t iloop;
- const t_ilist *il;
- int n, nmol;
+ gmx_mtop_ilistloop_t iloop;
+ const InteractionList *il;
+ int n, nmol;
n = 0;
iloop = gmx_mtop_ilistloop_init(mtop);
while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
{
- n += nmol*il[ftype].nr/(1+NRAL(ftype));
+ n += nmol*il[ftype].size()/(1+NRAL(ftype));
}
if (mtop->bIntermolecularInteractions)
{
- n += mtop->intermolecular_ilist[ftype].nr/(1+NRAL(ftype));
+ n += (*mtop->intermolecular_ilist)[ftype].size()/(1+NRAL(ftype));
}
return n;
dest->index[dest->nr] = dest->nra;
}
-static void ilistcat(int ftype, t_ilist *dest, const t_ilist *src, int copies,
- int dnum, int snum)
+static void ilistcat(int ftype,
+ t_ilist *dest,
+ const InteractionList &src,
+ int copies,
+ int dnum,
+ int snum)
{
int nral, c, i, a;
nral = NRAL(ftype);
- dest->nalloc = dest->nr + copies*src->nr;
+ dest->nalloc = dest->nr + copies*src.size();
srenew(dest->iatoms, dest->nalloc);
for (c = 0; c < copies; c++)
{
- for (i = 0; i < src->nr; )
+ for (i = 0; i < src.size(); )
{
- dest->iatoms[dest->nr++] = src->iatoms[i++];
+ dest->iatoms[dest->nr++] = src.iatoms[i++];
for (a = 0; a < nral; a++)
{
- dest->iatoms[dest->nr++] = dnum + src->iatoms[i++];
+ dest->iatoms[dest->nr++] = dnum + src.iatoms[i++];
}
}
dnum += snum;
for (ftype = 0; ftype < F_NRE; ftype++)
{
if (bMergeConstr &&
- ftype == F_CONSTR && molt.ilist[F_CONSTRNC].nr > 0)
+ ftype == F_CONSTR && molt.ilist[F_CONSTRNC].size() > 0)
{
/* Merge all constrains into one ilist.
* This simplifies the constraint code.
*/
for (int mol = 0; mol < molb.nmol; mol++)
{
- ilistcat(ftype, &idef->il[F_CONSTR], &molt.ilist[F_CONSTR],
+ ilistcat(ftype, &idef->il[F_CONSTR], molt.ilist[F_CONSTR],
1, destnr + mol*srcnr, srcnr);
- ilistcat(ftype, &idef->il[F_CONSTR], &molt.ilist[F_CONSTRNC],
+ ilistcat(ftype, &idef->il[F_CONSTR], molt.ilist[F_CONSTRNC],
1, destnr + mol*srcnr, srcnr);
}
}
else if (!(bMergeConstr && ftype == F_CONSTRNC))
{
- ilistcat(ftype, &idef->il[ftype], &molt.ilist[ftype],
+ ilistcat(ftype, &idef->il[ftype], molt.ilist[ftype],
molb.nmol, destnr, srcnr);
}
}
{
for (ftype = 0; ftype < F_NRE; ftype++)
{
- ilistcat(ftype, &idef->il[ftype], &mtop->intermolecular_ilist[ftype],
+ ilistcat(ftype, &idef->il[ftype], (*mtop->intermolecular_ilist)[ftype],
1, 0, mtop->natoms);
}
}
struct t_atom;
struct t_atoms;
struct t_block;
-struct t_ilist;
+struct InteractionList;
struct t_symtab;
// TODO All of the functions taking a const gmx_mtop * are deprecated
* When at the end, destroys iloop and returns FALSE.
*/
gmx_bool
-gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop,
- const t_ilist **ilist_mol, int *nmol);
-
+gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop,
+ const InteractionList **ilist_mol,
+ int *nmol);
/* Abstract type for ilist loop over all ilists of all molecules */
typedef struct gmx_mtop_ilistloop_all *gmx_mtop_ilistloop_all_t;
* When at the end, destroys iloop and returns FALSE.
*/
gmx_bool
-gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop,
- const t_ilist **ilist_mol, int *atnr_offset);
+gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop,
+ const InteractionList **ilist_mol,
+ int *atnr_offset);
/* Returns the total number of interactions in the system of type ftype */
excls()
{
init_t_atoms(&atoms, 0, FALSE);
-
- for (int ftype = 0; ftype < F_NRE; ftype++)
- {
- ilist[ftype].nr = 0;
- ilist[ftype].nr_nonperturbed = 0;
- ilist[ftype].iatoms = nullptr;
- ilist[ftype].nalloc = 0;
- }
}
gmx_moltype_t::~gmx_moltype_t()
done_atom(&atoms);
done_block(&cgs);
done_blocka(&excls);
-
- for (int f = 0; f < F_NRE; f++)
- {
- sfree(ilist[f].iatoms);
- ilist[f].nalloc = 0;
- }
}
void done_gmx_groups_t(gmx_groups_t *g)
{
pr_ilist(fp, indent, interaction_function[j].longname,
mtop->ffparams.functype,
- &mtop->intermolecular_ilist[j],
+ &(*mtop->intermolecular_ilist)[j],
bShowNumbers, bShowParameters, mtop->ffparams.iparams);
}
}
char **name; /**< Name of the molecule type */
t_atoms atoms; /**< The atoms in this molecule */
- t_ilist ilist[F_NRE]; /**< Interaction list with local indices */
+ InteractionList ilist[F_NRE]; /**< Interaction list with local indices */
t_block cgs; /**< The charge groups */
t_blocka excls; /**< The exclusions */
+
+ /* TODO: Change ilist to InteractionLists */
};
/*! \brief Block of molecules of the same type, used in gmx_mtop_t */
~gmx_mtop_t();
- char **name; /* Name of the topology */
- gmx_ffparams_t ffparams;
- std::vector<gmx_moltype_t> moltype;
- std::vector<gmx_molblock_t> molblock;
- gmx_bool bIntermolecularInteractions; /* Are there intermolecular
- * interactions? */
- t_ilist *intermolecular_ilist; /* List of intermolecular interactions
- * using system wide atom indices,
- * either NULL or size F_NRE */
+ char **name; /* Name of the topology */
+ gmx_ffparams_t ffparams;
+ std::vector<gmx_moltype_t> moltype;
+ std::vector<gmx_molblock_t> molblock;
+ gmx_bool bIntermolecularInteractions; /* Are there intermolecular
+ * interactions? */
+ std::unique_ptr<InteractionLists> intermolecular_ilist; /* List of intermolecular interactions
+ * using system wide atom indices,
+ * either NULL or size F_NRE */
int natoms;
- int maxres_renum; /* Parameter for residue numbering */
- int maxresnr; /* The maximum residue number in moltype */
- t_atomtypes atomtypes; /* Atomtype properties */
- gmx_groups_t groups; /* Groups of atoms for different purposes */
- t_symtab symtab; /* The symbol table */
- bool haveMoleculeIndices; /* Tells whether we have valid molecule indices */
+ int maxres_renum; /* Parameter for residue numbering */
+ int maxresnr; /* The maximum residue number in moltype */
+ t_atomtypes atomtypes; /* Atomtype properties */
+ gmx_groups_t groups; /* Groups of atoms for different purposes */
+ t_symtab symtab; /* The symbol table */
+ bool haveMoleculeIndices; /* Tells whether we have valid molecule indices */
/* Derived data */
std::vector<MoleculeBlockIndices> moleculeBlockIndices; /* Indices for each molblock entry for fast lookup of atom properties */
#include "gromacs/topology/ifunc.h"
#include "gromacs/topology/topology.h"
+#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
/* Check perturbed charges for 1-4 interactions */
for (const gmx_molblock_t &molb : mtop->molblock)
{
- const t_atom *atom = mtop->moltype[molb.type].atoms.atom;
- const t_ilist &il = mtop->moltype[molb.type].ilist[F_LJ14];
- const int *ia = il.iatoms;
- for (int i = 0; i < il.nr; i += 3)
+ const t_atom *atom = mtop->moltype[molb.type].atoms.atom;
+ const InteractionList &il = mtop->moltype[molb.type].ilist[F_LJ14];
+ gmx::ArrayRef<const int> ia = il.iatoms;
+ for (int i = 0; i < il.size(); i += 3)
{
if (atom[ia[i+1]].q != atom[ia[i+1]].qB ||
atom[ia[i+2]].q != atom[ia[i+2]].qB)
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff