This replaces t_idef in gmx_localtop_t.
Note that t_idef is still used in t_topology.
Extends InteractList with functionality to add entries
in a simple and clean way.
Change-Id: Ib59012332b6ba09df1de5dc696438c9aba4e2f2d
/*! \brief Sort ltop->ilist when we are doing free energy. */
void dd_sort_local_top(gmx_domdec_t* dd, const t_mdatoms* mdatoms, gmx_localtop_t* ltop);
-/*! \brief Initialize local topology
- *
- * \param[in] top_global Reference to global topology.
- * \param[in,out] top Pointer to new local topology
- */
-void dd_init_local_top(const gmx_mtop_t& top_global, gmx_localtop_t* top);
-
/*! \brief Construct local state */
void dd_init_local_state(struct gmx_domdec_t* dd, const t_state* state_global, t_state* local_state);
std::unique_ptr<gmx::HashedMap<int>> ga2la; /**< Global to local communicated constraint atom only index */
/* Multi-threading stuff */
- int nthread; /**< Number of threads used for DD constraint setup */
- std::vector<t_ilist> ils; /**< Constraint ilist working arrays, size \p nthread */
+ int nthread; /**< Number of threads used for DD constraint setup */
+ std::vector<InteractionList> ils; /**< Constraint ilist working arrays, size \p nthread */
/* Buffers for requesting atoms */
std::vector<std::vector<int>> requestedGlobalAtomIndices; /**< Buffers for requesting global atom indices, one per thread */
gmx_bool bHomeConnect,
gmx_domdec_constraints_t* dc,
gmx_domdec_specat_comm_t* dcc,
- t_ilist* il_local,
+ InteractionList* il_local,
std::vector<int>* ireq)
{
if (!dc->gc_req[con_offset + con])
/* Add this non-home constraint to the list */
dc->con_gl.push_back(con_offset + con);
dc->con_nlocat.push_back(bHomeConnect ? 1 : 0);
- dc->gc_req[con_offset + con] = true;
- if (il_local->nr + 3 > il_local->nalloc)
- {
- il_local->nalloc = over_alloc_dd(il_local->nr + 3);
- srenew(il_local->iatoms, il_local->nalloc);
- }
- const int* iap = constr_iatomptr(ia1, ia2, con);
- il_local->iatoms[il_local->nr++] = iap[0];
- const int a1_gl = offset + iap[1];
- const int a2_gl = offset + iap[2];
+ dc->gc_req[con_offset + con] = true;
+ const int* iap = constr_iatomptr(ia1, ia2, con);
+ const int parameterType = iap[0];
+ const int a1_gl = offset + iap[1];
+ const int a2_gl = offset + iap[2];
+ std::array<int, 2> atoms;
/* The following indexing code can probably be optizimed */
if (const int* a_loc = ga2la.findHome(a1_gl))
{
- il_local->iatoms[il_local->nr++] = *a_loc;
+ atoms[0] = *a_loc;
}
else
{
/* We set this index later */
- il_local->iatoms[il_local->nr++] = -a1_gl - 1;
+ atoms[0] = -a1_gl - 1;
}
if (const int* a_loc = ga2la.findHome(a2_gl))
{
- il_local->iatoms[il_local->nr++] = *a_loc;
+ atoms[1] = *a_loc;
}
else
{
/* We set this index later */
- il_local->iatoms[il_local->nr++] = -a2_gl - 1;
+ atoms[1] = -a2_gl - 1;
}
+ il_local->push_back(parameterType, atoms);
dc->ncon++;
}
/* Check to not ask for the same atom more than once */
gmx::ArrayRef<const std::vector<int>> at2settle_mt,
int cg_start,
int cg_end,
- t_ilist* ils_local,
+ InteractionList* ils_local,
std::vector<int>* ireq)
{
const gmx_ga2la_t& ga2la = *dd->ga2la;
if (bAssign)
{
- if (ils_local->nr + 1 + nral > ils_local->nalloc)
- {
- ils_local->nalloc = over_alloc_dd(ils_local->nr + 1 + nral);
- srenew(ils_local->iatoms, ils_local->nalloc);
- }
-
- ils_local->iatoms[ils_local->nr++] = ia1[settle * 4];
-
+ const int parameterType = ia1[settle * 4];
+ std::array<int, 3> atoms;
for (int sa = 0; sa < nral; sa++)
{
if (const int* a_loc = ga2la.findHome(a_gls[sa]))
{
- ils_local->iatoms[ils_local->nr++] = *a_loc;
+ atoms[sa] = *a_loc;
}
else
{
- ils_local->iatoms[ils_local->nr++] = -a_gls[sa] - 1;
+ atoms[sa] = -a_gls[sa] - 1;
/* Add this non-home atom to the list */
ireq->push_back(a_gls[sa]);
/* A check on double atom requests is
*/
}
}
+ ils_local->push_back(parameterType, atoms);
}
}
}
const int* cginfo,
gmx::ArrayRef<const ListOfLists<int>> at2con_mt,
int nrec,
- t_ilist* ilc_local,
+ InteractionList* ilc_local,
std::vector<int>* ireq)
{
gmx_domdec_constraints_t* dc = dd->constraints;
{
dc->con_gl.push_back(con_offset + con);
dc->con_nlocat.push_back(2);
- if (ilc_local->nr + 3 > ilc_local->nalloc)
- {
- ilc_local->nalloc = over_alloc_dd(ilc_local->nr + 3);
- srenew(ilc_local->iatoms, ilc_local->nalloc);
- }
- const int b_lo = *a_loc;
- ilc_local->iatoms[ilc_local->nr++] = iap[0];
- ilc_local->iatoms[ilc_local->nr++] = (a_gl == iap[1] ? a : b_lo);
- ilc_local->iatoms[ilc_local->nr++] = (a_gl == iap[1] ? b_lo : a);
+ const int b_lo = *a_loc;
+ const int parameterType = iap[0];
+ std::array<int, 2> atoms;
+ atoms[0] = (a_gl == iap[1] ? a : b_lo);
+ atoms[1] = (a_gl == iap[1] ? b_lo : a);
+ ilc_local->push_back(parameterType, atoms);
dc->ncon++;
nhome++;
}
}
}
-int dd_make_local_constraints(gmx_domdec_t* dd,
- int at_start,
- const struct gmx_mtop_t* mtop,
- const int* cginfo,
- gmx::Constraints* constr,
- int nrec,
- t_ilist* il_local)
+int dd_make_local_constraints(gmx_domdec_t* dd,
+ int at_start,
+ const struct gmx_mtop_t* mtop,
+ const int* cginfo,
+ gmx::Constraints* constr,
+ int nrec,
+ gmx::ArrayRef<InteractionList> il_local)
{
gmx_domdec_constraints_t* dc;
- t_ilist * ilc_local, *ils_local;
+ InteractionList * ilc_local, *ils_local;
gmx::HashedMap<int>* ga2la_specat;
int at_end, i, j;
- t_iatom* iap;
// This code should not be called unless this condition is true,
// because that's the only time init_domdec_constraints is
ilc_local = &il_local[F_CONSTR];
ils_local = &il_local[F_SETTLE];
- dc->ncon = 0;
- ilc_local->nr = 0;
+ dc->ncon = 0;
gmx::ArrayRef<const ListOfLists<int>> at2con_mt;
std::vector<int>* ireq = nullptr;
+ ilc_local->clear();
if (dd->constraint_comm)
{
// TODO Perhaps gmx_domdec_constraints_t should keep a valid constr?
{
// TODO Perhaps gmx_domdec_constraints_t should keep a valid constr?
GMX_RELEASE_ASSERT(constr != nullptr, "Must have valid constraints object");
- at2settle_mt = constr->atom2settle_moltype();
- ils_local->nr = 0;
+ at2settle_mt = constr->atom2settle_moltype();
+ ils_local->clear();
}
if (at2settle_mt.empty())
if (thread >= t0_set)
{
- int cg0, cg1;
- t_ilist* ilst;
+ int cg0, cg1;
+ InteractionList* ilst;
/* Distribute the settle check+assignments over
* dc->nthread or dc->nthread-1 threads.
{
ilst = &dc->ils[thread];
}
- ilst->nr = 0;
+ ilst->clear();
std::vector<int>& ireqt = dc->requestedGlobalAtomIndices[thread];
if (thread > 0)
/* Combine the generate settles and requested indices */
for (int thread = 1; thread < dc->nthread; thread++)
{
- t_ilist* ilst;
- int ia;
-
if (thread > t0_set)
{
- ilst = &dc->ils[thread];
- if (ils_local->nr + ilst->nr > ils_local->nalloc)
- {
- ils_local->nalloc = over_alloc_large(ils_local->nr + ilst->nr);
- srenew(ils_local->iatoms, ils_local->nalloc);
- }
- for (ia = 0; ia < ilst->nr; ia++)
- {
- ils_local->iatoms[ils_local->nr + ia] = ilst->iatoms[ia];
- }
- ils_local->nr += ilst->nr;
+ ils_local->append(dc->ils[thread]);
}
const std::vector<int>& ireqt = dc->requestedGlobalAtomIndices[thread];
if (debug)
{
- fprintf(debug, "Settles: total %3d\n", ils_local->nr / 4);
+ fprintf(debug, "Settles: total %3d\n", ils_local->size() / 4);
}
}
ga2la_specat = dd->constraints->ga2la.get();
nral1 = 1 + NRAL(F_CONSTR);
- for (i = 0; i < ilc_local->nr; i += nral1)
+ for (i = 0; i < ilc_local->size(); i += nral1)
{
- iap = ilc_local->iatoms + i;
+ int* iap = ilc_local->iatoms.data() + i;
for (j = 1; j < nral1; j++)
{
if (iap[j] < 0)
}
nral1 = 1 + NRAL(F_SETTLE);
- for (i = 0; i < ils_local->nr; i += nral1)
+ for (i = 0; i < ils_local->size(); i += nral1)
{
- iap = ils_local->iatoms + i;
+ int* iap = ils_local->iatoms.data() + i;
for (j = 1; j < nral1; j++)
{
if (iap[j] < 0)
#ifndef GMX_DOMDEC_DOMDEC_CONSTRAINTS_H
#define GMX_DOMDEC_DOMDEC_CONSTRAINTS_H
+#include "gromacs/utility/arrayref.h"
+
namespace gmx
{
class Constraints;
struct gmx_domdec_t;
struct gmx_mtop_t;
-struct t_ilist;
+struct InteractionList;
/*! \brief Clears the local indices for the constraint communication setup */
void dd_clear_local_constraint_indices(gmx_domdec_t* dd);
/*! \brief Sets up communication and atom indices for all local+connected constraints */
-int dd_make_local_constraints(struct gmx_domdec_t* dd,
- int at_start,
- const struct gmx_mtop_t* mtop,
- const int* cginfo,
- gmx::Constraints* constr,
- int nrec,
- struct t_ilist* il_local);
+int dd_make_local_constraints(struct gmx_domdec_t* dd,
+ int at_start,
+ const struct gmx_mtop_t* mtop,
+ const int* cginfo,
+ gmx::Constraints* constr,
+ int nrec,
+ gmx::ArrayRef<InteractionList> il_local);
/*! \brief Initializes the data structures for constraint communication */
void init_domdec_constraints(gmx_domdec_t* dd, const gmx_mtop_t* mtop);
/*! \brief Struct for thread local work data for local topology generation */
struct thread_work_t
{
- t_idef idef; /**< Partial local topology */
+ /*! \brief Constructor
+ *
+ * \param[in] ffparams The interaction parameters, the lifetime of the created object should not exceed the lifetime of the passed parameters
+ */
+ thread_work_t(const gmx_ffparams_t& ffparams) : idef(ffparams) {}
+
+ InteractionDefinitions idef; /**< Partial local topology */
std::unique_ptr<VsitePbc> vsitePbc; /**< vsite PBC structure */
int nbonded; /**< The number of bondeds in this struct */
ListOfLists<int> excl; /**< List of exclusions */
*
* \note This function needs to be called on all ranks (contains a global summation)
*/
-static std::string print_missing_interactions_mb(t_commrec* cr,
- const gmx_reverse_top_t* rt,
- const char* moltypename,
- const reverse_ilist_t* ril,
- int a_start,
- int a_end,
- int nat_mol,
- int nmol,
- const t_idef* idef)
+static std::string print_missing_interactions_mb(t_commrec* cr,
+ const gmx_reverse_top_t* rt,
+ const char* moltypename,
+ const reverse_ilist_t* ril,
+ int a_start,
+ int a_end,
+ int nat_mol,
+ int nmol,
+ const InteractionDefinitions* idef)
{
int* assigned;
int nril_mol = ril->index[nat_mol];
{
if (dd_check_ftype(ftype, rt->bBCheck, rt->bConstr, rt->bSettle))
{
- int nral = NRAL(ftype);
- const t_ilist* il = &idef->il[ftype];
- const t_iatom* ia = il->iatoms;
- for (int i = 0; i < il->nr; i += 1 + nral)
+ int nral = NRAL(ftype);
+ const InteractionList* il = &idef->il[ftype];
+ const int* ia = il->iatoms.data();
+ for (int i = 0; i < il->size(); i += 1 + nral)
{
int a0 = gatindex[ia[1]];
/* Check if this interaction is in
}
/*! \brief Help print error output when interactions are missing */
-static void print_missing_interactions_atoms(const gmx::MDLogger& mdlog,
- t_commrec* cr,
- const gmx_mtop_t* mtop,
- const t_idef* idef)
+static void print_missing_interactions_atoms(const gmx::MDLogger& mdlog,
+ t_commrec* cr,
+ const gmx_mtop_t* mtop,
+ const InteractionDefinitions* idef)
{
const gmx_reverse_top_t* rt = cr->dd->reverse_top;
for (ftype = 0; ftype < F_NRE; ftype++)
{
nral = NRAL(ftype);
- cl[ftype] = top_local->idef.il[ftype].nr / (1 + nral);
+ cl[ftype] = top_local->idef.il[ftype].size() / (1 + nral);
}
gmx_sumi(F_NRE, cl, cr);
rt.mbi.push_back(mbi);
}
- rt.th_work.resize(gmx_omp_nthreads_get(emntDomdec));
+ for (int th = 0; th < gmx_omp_nthreads_get(emntDomdec); th++)
+ {
+ rt.th_work.emplace_back(mtop->ffparams);
+ }
return rt;
}
int a_gl,
int a_mol,
const t_iatom* iatoms,
- t_ilist* il)
+ InteractionList* il)
{
- t_iatom* liatoms;
-
- if (il->nr + 1 + nral > il->nalloc)
- {
- il->nalloc = over_alloc_large(il->nr + 1 + nral);
- srenew(il->iatoms, il->nalloc);
- }
- liatoms = il->iatoms + il->nr;
- il->nr += 1 + nral;
-
/* Copy the type */
tiatoms[0] = iatoms[0];
// Note that ga2la_get_home always sets the third parameter if
// it returns TRUE
}
- for (int k = 0; k < 1 + nral; k++)
- {
- liatoms[k] = tiatoms[k];
- }
-}
-
-/*! \brief Store a bonded interaction at the end of \p il */
-static inline void add_ifunc(int nral, const t_iatom* tiatoms, t_ilist* il)
-{
- t_iatom* liatoms;
- int k;
-
- if (il->nr + 1 + nral > il->nalloc)
- {
- il->nalloc = over_alloc_large(il->nr + 1 + nral);
- srenew(il->iatoms, il->nalloc);
- }
- liatoms = il->iatoms + il->nr;
- for (k = 0; k <= nral; k++)
- {
- liatoms[k] = tiatoms[k];
- }
- il->nr += 1 + nral;
+ il->push_back(tiatoms[0], nral, tiatoms + 1);
}
/*! \brief Store a position restraint in idef and iatoms, complex because the parameters are different for each entry */
-static void add_posres(int mol,
- int a_mol,
- int numAtomsInMolecule,
- const gmx_molblock_t* molb,
- t_iatom* iatoms,
- const t_iparams* ip_in,
- t_idef* idef)
+static void add_posres(int mol,
+ int a_mol,
+ int numAtomsInMolecule,
+ const gmx_molblock_t* molb,
+ t_iatom* iatoms,
+ const t_iparams* ip_in,
+ InteractionDefinitions* idef)
{
- int n, a_molb;
- t_iparams* ip;
-
/* This position restraint has not been added yet,
* so it's index is the current number of position restraints.
*/
- n = idef->il[F_POSRES].nr / 2;
- if (n + 1 > idef->iparams_posres_nalloc)
- {
- idef->iparams_posres_nalloc = over_alloc_dd(n + 1);
- srenew(idef->iparams_posres, idef->iparams_posres_nalloc);
- }
- ip = &idef->iparams_posres[n];
- /* Copy the force constants */
- *ip = ip_in[iatoms[0]];
+ const int n = idef->il[F_POSRES].size() / 2;
/* Get the position restraint coordinates from the molblock */
- a_molb = mol * numAtomsInMolecule + a_mol;
+ const int a_molb = mol * numAtomsInMolecule + a_mol;
GMX_ASSERT(a_molb < ssize(molb->posres_xA),
"We need a sufficient number of position restraint coordinates");
- ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
- ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
- ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
+ /* Copy the force constants */
+ t_iparams ip = ip_in[iatoms[0]];
+ ip.posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
+ ip.posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
+ ip.posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
if (!molb->posres_xB.empty())
{
- ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
- ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
- ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
+ ip.posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
+ ip.posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
+ ip.posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
}
else
{
- ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
- ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
- ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
+ ip.posres.pos0B[XX] = ip.posres.pos0A[XX];
+ ip.posres.pos0B[YY] = ip.posres.pos0A[YY];
+ ip.posres.pos0B[ZZ] = ip.posres.pos0A[ZZ];
}
- /* Set the parameter index for idef->iparams_posre */
+ /* Set the parameter index for idef->iparams_posres */
iatoms[0] = n;
+ idef->iparams_posres.push_back(ip);
+
+ GMX_ASSERT(int(idef->iparams_posres.size()) == n + 1,
+ "The index of the parameter type should match n");
}
/*! \brief Store a flat-bottomed position restraint in idef and iatoms, complex because the parameters are different for each entry */
-static void add_fbposres(int mol,
- int a_mol,
- int numAtomsInMolecule,
- const gmx_molblock_t* molb,
- t_iatom* iatoms,
- const t_iparams* ip_in,
- t_idef* idef)
+static void add_fbposres(int mol,
+ int a_mol,
+ int numAtomsInMolecule,
+ const gmx_molblock_t* molb,
+ t_iatom* iatoms,
+ const t_iparams* ip_in,
+ InteractionDefinitions* idef)
{
- int n, a_molb;
- t_iparams* ip;
-
/* This flat-bottom position restraint has not been added yet,
* so it's index is the current number of position restraints.
*/
- n = idef->il[F_FBPOSRES].nr / 2;
- if (n + 1 > idef->iparams_fbposres_nalloc)
- {
- idef->iparams_fbposres_nalloc = over_alloc_dd(n + 1);
- srenew(idef->iparams_fbposres, idef->iparams_fbposres_nalloc);
- }
- ip = &idef->iparams_fbposres[n];
- /* Copy the force constants */
- *ip = ip_in[iatoms[0]];
+ const int n = idef->il[F_FBPOSRES].size() / 2;
/* Get the position restraint coordinats from the molblock */
- a_molb = mol * numAtomsInMolecule + a_mol;
+ const int a_molb = mol * numAtomsInMolecule + a_mol;
GMX_ASSERT(a_molb < ssize(molb->posres_xA),
"We need a sufficient number of position restraint coordinates");
+ /* Copy the force constants */
+ t_iparams ip = ip_in[iatoms[0]];
/* Take reference positions from A position of normal posres */
- ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
- ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
- ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
+ ip.fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
+ ip.fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
+ ip.fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
/* Note: no B-type for flat-bottom posres */
- /* Set the parameter index for idef->iparams_posre */
+ /* Set the parameter index for idef->iparams_fbposres */
iatoms[0] = n;
+ idef->iparams_fbposres.push_back(ip);
+
+ GMX_ASSERT(int(idef->iparams_fbposres.size()) == n + 1,
+ "The index of the parameter type should match n");
}
/*! \brief Store a virtual site interaction, complex because of PBC and recursion */
int a_gl,
int a_mol,
const t_iatom* iatoms,
- t_idef* idef)
+ InteractionDefinitions* idef)
{
int k;
t_iatom tiatoms[1 + MAXATOMLIST];
}
/*! \brief Append t_idef structures 1 to nsrc in src to *dest */
-static void combine_idef(t_idef* dest, gmx::ArrayRef<const thread_work_t> src)
+static void combine_idef(InteractionDefinitions* dest, gmx::ArrayRef<const thread_work_t> src)
{
int ftype;
int n = 0;
for (gmx::index s = 1; s < src.ssize(); s++)
{
- n += src[s].idef.il[ftype].nr;
+ n += src[s].idef.il[ftype].size();
}
if (n > 0)
{
- t_ilist* ild;
-
- ild = &dest->il[ftype];
-
- if (ild->nr + n > ild->nalloc)
- {
- ild->nalloc = over_alloc_large(ild->nr + n);
- srenew(ild->iatoms, ild->nalloc);
- }
-
for (gmx::index s = 1; s < src.ssize(); s++)
{
- const t_ilist& ils = src[s].idef.il[ftype];
-
- for (int i = 0; i < ils.nr; i++)
- {
- ild->iatoms[ild->nr + i] = ils.iatoms[i];
- }
-
- ild->nr += ils.nr;
+ dest->il[ftype].append(src[s].idef.il[ftype]);
}
/* Position restraints need an additional treatment */
if (ftype == F_POSRES || ftype == F_FBPOSRES)
{
- int nposres = dest->il[ftype].nr / 2;
- // TODO: Simplify this code using std::vector
- t_iparams*& iparams_dest =
+ int nposres = dest->il[ftype].size() / 2;
+ std::vector<t_iparams>& iparams_dest =
(ftype == F_POSRES ? dest->iparams_posres : dest->iparams_fbposres);
- int& posres_nalloc = (ftype == F_POSRES ? dest->iparams_posres_nalloc
- : dest->iparams_fbposres_nalloc);
- if (nposres > posres_nalloc)
- {
- posres_nalloc = over_alloc_large(nposres);
- srenew(iparams_dest, posres_nalloc);
- }
/* Set nposres to the number of original position restraints in dest */
for (gmx::index s = 1; s < src.ssize(); s++)
{
- nposres -= src[s].idef.il[ftype].nr / 2;
+ nposres -= src[s].idef.il[ftype].size() / 2;
}
for (gmx::index s = 1; s < src.ssize(); s++)
{
- const t_iparams* iparams_src = (ftype == F_POSRES ? src[s].idef.iparams_posres
- : src[s].idef.iparams_fbposres);
+ const std::vector<t_iparams>& iparams_src =
+ (ftype == F_POSRES ? src[s].idef.iparams_posres : src[s].idef.iparams_fbposres);
+ iparams_dest.insert(iparams_dest.end(), iparams_src.begin(), iparams_src.end());
- for (int i = 0; i < src[s].idef.il[ftype].nr / 2; i++)
+ /* Correct the indices into iparams_posres */
+ for (int i = 0; i < src[s].idef.il[ftype].size() / 2; i++)
{
/* Correct the index into iparams_posres */
dest->il[ftype].iatoms[nposres * 2] = nposres;
- /* Copy the position restraint force parameters */
- iparams_dest[nposres] = iparams_src[i];
nposres++;
}
}
+ GMX_RELEASE_ASSERT(
+ int(iparams_dest.size()) == nposres,
+ "The number of parameters should match the number of restraints");
}
}
}
t_pbc* pbc_null,
rvec* cg_cm,
const t_iparams* ip_in,
- t_idef* idef,
+ InteractionDefinitions* idef,
int iz,
gmx_bool bBCheck,
int* nbonded_local)
if (bUse)
{
/* Add this interaction to the local topology */
- add_ifunc(nral, tiatoms, &idef->il[ftype]);
+ idef->il[ftype].push_back(tiatoms[0], nral, tiatoms + 1);
/* Sum so we can check in global_stat
* if we have everything.
*/
t_pbc* pbc_null,
rvec* cg_cm,
const t_iparams* ip_in,
- t_idef* idef,
+ InteractionDefinitions* idef,
int izone,
const gmx::Range<int>& atomRange)
{
"The number of exclusion list should match the number of atoms in the range");
}
-/*! \brief Clear a t_idef struct */
-static void clear_idef(t_idef* idef)
-{
- int ftype;
-
- /* Clear the counts */
- for (ftype = 0; ftype < F_NRE; ftype++)
- {
- idef->il[ftype].nr = 0;
- }
-}
-
/*! \brief Generate and store all required local bonded interactions in \p idef and local exclusions in \p lexcls */
-static int make_local_bondeds_excls(gmx_domdec_t* dd,
- gmx_domdec_zones_t* zones,
- const gmx_mtop_t* mtop,
- const int* cginfo,
- gmx_bool bRCheckMB,
- ivec rcheck,
- gmx_bool bRCheck2B,
- real rc,
- t_pbc* pbc_null,
- rvec* cg_cm,
- t_idef* idef,
- ListOfLists<int>* lexcls,
- int* excl_count)
+static int make_local_bondeds_excls(gmx_domdec_t* dd,
+ gmx_domdec_zones_t* zones,
+ const gmx_mtop_t* mtop,
+ const int* cginfo,
+ gmx_bool bRCheckMB,
+ ivec rcheck,
+ gmx_bool bRCheck2B,
+ real rc,
+ t_pbc* pbc_null,
+ rvec* cg_cm,
+ InteractionDefinitions* idef,
+ ListOfLists<int>* lexcls,
+ int* excl_count)
{
int nzone_bondeds;
int cg0, cg1;
rc2 = rc * rc;
/* Clear the counts */
- clear_idef(idef);
+ idef->clear();
nbonded_local = 0;
lexcls->clear();
{
try
{
- int cg0t, cg1t;
- t_idef* idef_t;
+ int cg0t, cg1t;
+ InteractionDefinitions* idef_t;
cg0t = cg0 + ((cg1 - cg0) * thread) / numThreads;
cg1t = cg0 + ((cg1 - cg0) * (thread + 1)) / numThreads;
else
{
idef_t = &rt->th_work[thread].idef;
- clear_idef(idef_t);
+ idef_t->clear();
}
rt->th_work[thread].nbonded = make_bondeds_zone(
dd, zones, mtop->molblock, bRCheckMB, rcheck, bRCheck2B, rc2, pbc_null,
- cg_cm, idef->iparams, idef_t, izone, gmx::Range<int>(cg0t, cg1t));
+ cg_cm, idef->iparams.data(), idef_t, izone, gmx::Range<int>(cg0t, cg1t));
if (izone < numIZonesForExclusions)
{
* we can only do this when we have the charge arrays.
*/
ltop->idef.ilsort = ilsortUNKNOWN;
-
- ltop->atomtypes = mtop.atomtypes;
}
void dd_sort_local_top(gmx_domdec_t* dd, const t_mdatoms* mdatoms, gmx_localtop_t* ltop)
}
}
-void dd_init_local_top(const gmx_mtop_t& top_global, gmx_localtop_t* top)
-{
- /* TODO: Get rid of the const casts below, e.g. by using a reference */
- top->idef.ntypes = top_global.ffparams.numTypes();
- top->idef.atnr = top_global.ffparams.atnr;
- top->idef.functype = const_cast<t_functype*>(top_global.ffparams.functype.data());
- top->idef.iparams = const_cast<t_iparams*>(top_global.ffparams.iparams.data());
- top->idef.fudgeQQ = top_global.ffparams.fudgeQQ;
- top->idef.cmap_grid = new gmx_cmap_t;
- *top->idef.cmap_grid = top_global.ffparams.cmap_grid;
-
- top->idef.ilsort = ilsortUNKNOWN;
- top->useInDomainDecomp_ = true;
-}
-
void dd_init_local_state(gmx_domdec_t* dd, const t_state* state_global, t_state* state_local)
{
int buf[NITEM_DD_INIT_LOCAL_STATE];
bool hasVsite = false;
for (int i = 0; i < F_NRE; i++)
{
- if ((interaction_function[i].flags & IF_VSITE) && molt.ilist[i].size() > 0)
+ if ((interaction_function[i].flags & IF_VSITE) && !molt.ilist[i].empty())
{
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.data(), ilist, PbcType::No,
+ construct_vsites(nullptr, xs, 0.0, nullptr, ffparams->iparams, molt->ilist, PbcType::No,
TRUE, nullptr, nullptr);
}
}
}
}
-int dd_make_local_vsites(gmx_domdec_t* dd, int at_start, t_ilist* lil)
+int dd_make_local_vsites(gmx_domdec_t* dd, int at_start, gmx::ArrayRef<InteractionList> lil)
{
std::vector<int>& ireq = dd->vsite_requestedGlobalAtomIndices;
gmx::HashedMap<int>* ga2la_specat = dd->ga2la_vsite;
{
if (interaction_function[ftype].flags & IF_VSITE)
{
- const int nral = NRAL(ftype);
- const t_ilist& lilf = lil[ftype];
- for (int i = 0; i < lilf.nr; i += 1 + nral)
+ const int nral = NRAL(ftype);
+ const InteractionList& lilf = lil[ftype];
+ for (int i = 0; i < lilf.size(); i += 1 + nral)
{
- const t_iatom* iatoms = lilf.iatoms + i;
+ const int* iatoms = lilf.iatoms.data() + i;
/* Check if we have the other atoms */
for (int j = 1; j < 1 + nral; j++)
{
{
if (interaction_function[ftype].flags & IF_VSITE)
{
- const int nral = NRAL(ftype);
- t_ilist& lilf = lil[ftype];
- for (int i = 0; i < lilf.nr; i += 1 + nral)
+ const int nral = NRAL(ftype);
+ InteractionList& lilf = lil[ftype];
+ for (int i = 0; i < lilf.size(); i += 1 + nral)
{
- t_iatom* iatoms = lilf.iatoms + i;
+ t_iatom* iatoms = lilf.iatoms.data() + i;
for (int j = 1; j < 1 + nral; j++)
{
if (iatoms[j] < 0)
#ifndef GMX_DOMDEC_DOMDEC_VSITE_H
#define GMX_DOMDEC_DOMDEC_VSITE_H
-struct t_ilist;
+#include "gromacs/utility/arrayref.h"
+
struct gmx_domdec_t;
+struct InteractionList;
/*! \brief Clears the local indices for the virtual site communication setup */
void dd_clear_local_vsite_indices(struct gmx_domdec_t* dd);
/*! \brief Sets up communication and atom indices for all local vsites */
-int dd_make_local_vsites(struct gmx_domdec_t* dd, int at_start, struct t_ilist* lil);
+int dd_make_local_vsites(struct gmx_domdec_t* dd, int at_start, gmx::ArrayRef<InteractionList> lil);
/*! \brief Initializes the data structures for virtual site communication */
void init_domdec_vsites(struct gmx_domdec_t* dd, int n_intercg_vsite);
{
GMX_ASSERT(graph != nullptr, "We use a graph with PBC (no periodic mols) and without DD");
- *graph = mk_graph(nullptr, &(top->idef), 0, top_global.natoms, FALSE, FALSE);
+ *graph = mk_graph(nullptr, top->idef, 0, top_global.natoms, FALSE, FALSE);
}
else if (graph != nullptr)
{
if (constr)
{
- constr->setConstraints(*top, *mdatoms);
+ constr->setConstraints(top, *mdatoms);
}
}
else
{
do_ilist(serializer, &ilist);
- if (file_version < 78 && j == F_SETTLE && ilist.size() > 0)
+ if (file_version < 78 && j == F_SETTLE && !ilist.empty())
{
add_settle_atoms(&ilist);
}
{
const InteractionList& il = (*ilist)[F_DISRES];
- if (il.size() > 0)
+ if (!il.empty())
{
gmx::ArrayRef<const int> a = il.iatoms;
int npair = 0;
{
if (tpx->fileVersion < 57)
{
- if (mtop->moltype[0].ilist[F_DISRES].size() > 0)
+ if (!mtop->moltype[0].ilist[F_DISRES].empty())
{
ir->eDisre = edrSimple;
}
fprintf(fp, "\n");
}
-static void check_viol(FILE* log,
- t_ilist* disres,
- t_iparams forceparams[],
- rvec x[],
- rvec4 f[],
- t_pbc* pbc,
- t_graph* g,
- t_dr_result dr[],
- int clust_id,
- int isize,
- const int index[],
- real vvindex[],
- t_fcdata* fcd)
+static void check_viol(FILE* log,
+ const InteractionList& disres,
+ gmx::ArrayRef<const t_iparams> forceparams,
+ rvec x[],
+ rvec4 f[],
+ t_pbc* pbc,
+ t_graph* g,
+ t_dr_result dr[],
+ int clust_id,
+ int isize,
+ const int index[],
+ real vvindex[],
+ t_fcdata* fcd)
{
- t_iatom* forceatoms;
int i, j, nat, n, type, nviol, ndr, label;
real rt, mviol, tviol, viol, lam, dvdl, drt;
rvec* fshift;
{
reset5();
}
- forceatoms = disres->iatoms;
+ gmx::ArrayRef<const int> forceatoms = disres.iatoms;
for (j = 0; (j < isize); j++)
{
vvindex[j] = 0;
// The label for a distance restraint should be at most one larger
// than the previous label.
int label_old = forceparams[forceatoms[0]].disres.label;
- for (i = 0; (i < disres->nr); i += nat)
+ for (i = 0; (i < disres.size()); i += nat)
{
type = forceatoms[i];
label = forceparams[type].disres.label;
}
// Get offset for label index
label_old = forceparams[forceatoms[0]].disres.label;
- for (i = 0; (i < disres->nr);)
+ for (i = 0; (i < disres.size());)
{
type = forceatoms[i];
n = 0;
do
{
n += nat;
- } while (((i + n) < disres->nr)
+ } while (((i + n) < disres.size())
&& (forceparams[forceatoms[i + n]].disres.label == label + label_old));
calc_disres_R_6(nullptr, nullptr, n, &forceatoms[i], x, pbc, fcd, nullptr);
dr[clust_id].aver_6[ndr] += fcd->disres.Rt_6[label];
snew(fshift, SHIFTS);
- ta_disres(n, &forceatoms[i], forceparams, x, f, fshift, pbc, g, lam, &dvdl, nullptr, fcd, nullptr);
+ ta_disres(n, &forceatoms[i], forceparams.data(), x, f, fshift, pbc, g, lam, &dvdl, nullptr,
+ fcd, nullptr);
sfree(fshift);
viol = fcd->disres.sumviol;
if (bFirst)
{
- fprintf(stderr, "\nThere are %d restraints and %d pairs\n", ndr, disres->nr / nat);
+ fprintf(stderr, "\nThere are %d restraints and %d pairs\n", ndr, disres.size() / nat);
bFirst = FALSE;
}
if (ntop)
return bIC;
}
-static void dump_stats(FILE* log,
- int nsteps,
- const t_disresdata& dd,
- const t_ilist* disres,
- t_iparams ip[],
- t_dr_result* dr,
- int isize,
- int index[],
- t_atoms* atoms)
+static void dump_stats(FILE* log,
+ int nsteps,
+ const t_disresdata& dd,
+ const InteractionList& disres,
+ gmx::ArrayRef<const t_iparams> ip,
+ t_dr_result* dr,
+ int isize,
+ int index[],
+ t_atoms* atoms)
{
t_dr_stats* drs;
fprintf(log, "++++++++++++++ STATISTICS ++++++++++++++++++++++++\n");
snew(drs, dd.nres);
const int nra = interaction_function[F_DISRES].nratoms + 1;
- for (int j = 0; j < disres->nr; j += nra)
+ for (int j = 0; j < disres.size(); j += nra)
{
// Note that the restraint i can be used by multiple pairs
- const int i = disres->iatoms[j] - dd.type_min;
+ const int i = disres.iatoms[j] - dd.type_min;
GMX_RELEASE_ASSERT(i >= 0 && i < dd.nres, "The restraint index should be in range");
- drs[i].label = ip[disres->iatoms[j]].disres.label;
+ drs[i].label = ip[disres.iatoms[j]].disres.label;
drs[i].bCore = is_core(drs[i].label, isize, index);
- drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
+ drs[i].up1 = ip[disres.iatoms[j]].disres.up1;
drs[i].r = dr->aver1[i] / nsteps;
drs[i].rT3 = gmx::invcbrt(dr->aver_3[i] / nsteps);
drs[i].rT6 = gmx::invsixthroot(dr->aver_6[i] / nsteps);
drs[i].violT6 = std::max(0.0, static_cast<double>(drs[i].rT6 - drs[i].up1));
if (atoms)
{
- int j1 = disres->iatoms[j + 1];
- int j2 = disres->iatoms[j + 2];
+ int j1 = disres.iatoms[j + 1];
+ int j2 = disres.iatoms[j + 2];
atoms->pdbinfo[j1].bfac += drs[i].violT3 * 5;
atoms->pdbinfo[j2].bfac += drs[i].violT3 * 5;
}
sfree(drs);
}
-static void dump_clust_stats(FILE* fp,
- const t_disresdata& dd,
- const t_ilist* disres,
- t_iparams ip[],
- t_blocka* clust,
- t_dr_result dr[],
- char* clust_name[],
- int isize,
- int index[])
+static void dump_clust_stats(FILE* fp,
+ const t_disresdata& dd,
+ const InteractionList& disres,
+ gmx::ArrayRef<const t_iparams> ip,
+ t_blocka* clust,
+ t_dr_result dr[],
+ char* clust_name[],
+ int isize,
+ int index[])
{
int k, nra, mmm = 0;
double sumV, maxV, sumVT3, sumVT6, maxVT3, maxVT6;
for (int j = 0; j < dd.nres; j += nra)
{
// Note that the restraint i can be used by multiple pairs
- const int i = disres->iatoms[j] - dd.type_min;
+ const int i = disres.iatoms[j] - dd.type_min;
if (restraintHasBeenProcessed[i])
{
continue;
}
- drs[i].label = ip[disres->iatoms[j]].disres.label;
+ drs[i].label = ip[disres.iatoms[j]].disres.label;
drs[i].bCore = is_core(drs[i].label, isize, index);
- drs[i].up1 = ip[disres->iatoms[j]].disres.up1;
+ drs[i].up1 = ip[disres.iatoms[j]].disres.up1;
drs[i].r = dr[k].aver1[i] / dr[k].nframes;
if ((dr[k].aver_3[i] <= 0) || !std::isfinite(dr[k].aver_3[i]))
{
dr->averv = 0;
}
-static void dump_disre_matrix(const char* fn,
- t_dr_result* dr,
- int ndr,
- int nsteps,
- t_idef* idef,
- const gmx_mtop_t* mtop,
- real max_dr,
- int nlevels,
- gmx_bool bThird)
+static void dump_disre_matrix(const char* fn,
+ t_dr_result* dr,
+ int ndr,
+ int nsteps,
+ const InteractionDefinitions& idef,
+ const gmx_mtop_t* mtop,
+ real max_dr,
+ int nlevels,
+ gmx_bool bThird)
{
FILE* fp;
int* resnr;
snew(matrix[i], n_res);
}
nratoms = interaction_function[F_DISRES].nratoms;
- nra = (idef->il[F_DISRES].nr / (nratoms + 1));
+ nra = (idef.il[F_DISRES].size() / (nratoms + 1));
snew(ptr, nra + 1);
index = 0;
nlabel = 0;
ptr[0] = 0;
snew(w_dr, ndr);
- for (i = 0; (i < idef->il[F_DISRES].nr); i += nratoms + 1)
+ for (i = 0; (i < idef.il[F_DISRES].size()); i += nratoms + 1)
{
- tp = idef->il[F_DISRES].iatoms[i];
- label = idef->iparams[tp].disres.label;
+ tp = idef.il[F_DISRES].iatoms[i];
+ label = idef.iparams[tp].disres.label;
if (label != index)
{
{
for (j = ptr[i]; (j < ptr[i + 1]); j += nratoms + 1)
{
- tp = idef->il[F_DISRES].iatoms[j];
- ai = idef->il[F_DISRES].iatoms[j + 1];
- aj = idef->il[F_DISRES].iatoms[j + 2];
+ tp = idef.il[F_DISRES].iatoms[j];
+ ai = idef.il[F_DISRES].iatoms[j + 1];
+ aj = idef.il[F_DISRES].iatoms[j + 2];
ri = resnr[ai];
rj = resnr[aj];
{
fprintf(debug, "DR %d, atoms %d, %d, distance %g\n", i, ai, aj, rav);
}
- rviol = std::max(0.0_real, rav - idef->iparams[tp].disres.up1);
+ rviol = std::max(0.0_real, rav - idef.iparams[tp].disres.up1);
matrix[ri][rj] += w_dr[i] * rviol;
matrix[rj][ri] += w_dr[i] * rviol;
hi = std::max(hi, matrix[ri][rj]);
"Use inverse third power averaging or linear for matrix output" }
};
- FILE *out = nullptr, *aver = nullptr, *numv = nullptr, *maxxv = nullptr, *xvg = nullptr;
- gmx_localtop_t top;
- t_fcdata fcd;
- t_graph* g;
- int i, j, kkk;
- t_trxstatus* status;
- real t;
- rvec * x, *xav = nullptr;
- rvec4* f;
- matrix box;
- gmx_bool bPDB;
- int isize;
- int * index = nullptr, *ind_fit = nullptr;
- char* grpname;
+ FILE * out = nullptr, *aver = nullptr, *numv = nullptr, *maxxv = nullptr, *xvg = nullptr;
+ t_fcdata fcd;
+ t_graph* g;
+ int i, j, kkk;
+ t_trxstatus* status;
+ real t;
+ rvec * x, *xav = nullptr;
+ rvec4* f;
+ matrix box;
+ gmx_bool bPDB;
+ int isize;
+ int * index = nullptr, *ind_fit = nullptr;
+ char* grpname;
t_cluster_ndx* clust = nullptr;
t_dr_result dr, *dr_clust = nullptr;
char** leg;
atoms->havePdbInfo = TRUE;
}
+ gmx_localtop_t top(topInfo.mtop()->ffparams);
gmx_mtop_generate_local_top(*topInfo.mtop(), &top, ir->efep != efepNO);
g = nullptr;
}
else
{
- g = mk_graph(fplog, &top.idef, 0, ntopatoms, FALSE, FALSE);
+ g = mk_graph(fplog, top.idef, 0, ntopatoms, FALSE, FALSE);
}
}
update_mdatoms(mdAtoms->mdatoms(), ir->fepvals->init_lambda);
if (ir->pbcType != PbcType::No)
{
- gpbc = gmx_rmpbc_init(&top.idef, ir->pbcType, natoms);
+ gpbc = gmx_rmpbc_init(top.idef, ir->pbcType, natoms);
}
j = 0;
}
my_clust = clust->inv_clust[j];
range_check(my_clust, 0, clust->clust->nr);
- check_viol(fplog, &(top.idef.il[F_DISRES]), top.idef.iparams, x, f, pbc_null, g,
- dr_clust, my_clust, isize, index, vvindex, &fcd);
+ check_viol(fplog, top.idef.il[F_DISRES], top.idef.iparams, x, f, pbc_null, g, dr_clust,
+ my_clust, isize, index, vvindex, &fcd);
}
else
{
- check_viol(fplog, &(top.idef.il[F_DISRES]), top.idef.iparams, x, f, pbc_null, g, &dr, 0,
+ check_viol(fplog, top.idef.il[F_DISRES], top.idef.iparams, x, f, pbc_null, g, &dr, 0,
isize, index, vvindex, &fcd);
}
if (bPDB)
if (clust)
{
- dump_clust_stats(fplog, fcd.disres, &(top.idef.il[F_DISRES]), top.idef.iparams,
- clust->clust, dr_clust, clust->grpname, isize, index);
+ dump_clust_stats(fplog, fcd.disres, top.idef.il[F_DISRES], top.idef.iparams, clust->clust,
+ dr_clust, clust->grpname, isize, index);
}
else
{
- dump_stats(fplog, j, fcd.disres, &(top.idef.il[F_DISRES]), top.idef.iparams, &dr, isize,
- index, bPDB ? atoms.get() : nullptr);
+ dump_stats(fplog, j, fcd.disres, top.idef.il[F_DISRES], top.idef.iparams, &dr, isize, index,
+ bPDB ? atoms.get() : nullptr);
if (bPDB)
{
write_sto_conf(opt2fn("-q", NFILE, fnm), "Coloured by average violation in Angstrom",
atoms.get(), xav, nullptr, ir->pbcType, box);
}
- dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres, j, &top.idef,
+ dump_disre_matrix(opt2fn_null("-x", NFILE, fnm), &dr, fcd.disres.nres, j, top.idef,
topInfo.mtop(), max_dr, nlevels, bThird);
xvgrclose(out);
xvgrclose(aver);
static int get_bounds(real** bounds, int** index, int** dr_pair, int* npairs, gmx_localtop_t* top)
{
- t_functype* functype;
- t_iparams* ip;
- int i, j, k, type, ftype, natom;
- t_ilist* disres;
- t_iatom* iatom;
- real* b;
- int * ind, *pair;
- int nb, label1;
-
- functype = top->idef.functype;
- ip = top->idef.iparams;
+ int i, j, k, type, ftype, natom;
+ real* b;
+ int * ind, *pair;
+ int nb, label1;
+
+ gmx::ArrayRef<const t_functype> functype = top->idef.functype;
+ gmx::ArrayRef<const t_iparams> iparams = top->idef.iparams;
/* Count how many distance restraint there are... */
- nb = top->idef.il[F_DISRES].nr;
+ nb = top->idef.il[F_DISRES].size();
if (nb == 0)
{
gmx_fatal(FARGS, "No distance restraints in topology!\n");
/* Fill the bound array */
nb = 0;
- for (i = 0; (i < top->idef.ntypes); i++)
+ for (gmx::index i = 0; i < functype.ssize(); i++)
{
ftype = functype[i];
if (ftype == F_DISRES)
{
- label1 = ip[i].disres.label;
- b[nb] = ip[i].disres.up1;
+ label1 = iparams[i].disres.label;
+ b[nb] = iparams[i].disres.up1;
ind[nb] = label1;
nb++;
}
*bounds = b;
/* Fill the index array */
- label1 = -1;
- disres = &(top->idef.il[F_DISRES]);
- iatom = disres->iatoms;
- for (i = j = k = 0; (i < disres->nr);)
+ label1 = -1;
+ const InteractionList& disres = top->idef.il[F_DISRES];
+ gmx::ArrayRef<const int> iatom = disres.iatoms;
+ for (i = j = k = 0; (i < disres.size());)
{
type = iatom[i];
- ftype = top->idef.functype[type];
+ ftype = functype[type];
natom = interaction_function[ftype].nratoms + 1;
- if (label1 != top->idef.iparams[type].disres.label)
+ if (label1 != iparams[type].disres.label)
{
pair[j] = k;
- label1 = top->idef.iparams[type].disres.label;
+ label1 = iparams[type].disres.label;
j++;
}
k++;
FILE /* *out = NULL,*/ *out_disre = nullptr, *fp_pairs = nullptr, *fort = nullptr,
*fodt = nullptr, *foten = nullptr;
- ener_file_t fp;
- int timecheck = 0;
- gmx_localtop_t top;
- gmx_enxnm_t* enm = nullptr;
- t_enxframe fr;
- int nre, teller, teller_disre;
- int nor = 0, nex = 0, norfr = 0, enx_i = 0;
+ ener_file_t fp;
+ int timecheck = 0;
+ gmx_enxnm_t* enm = nullptr;
+ t_enxframe fr;
+ int nre, teller, teller_disre;
+ int nor = 0, nex = 0, norfr = 0, enx_i = 0;
real *bounds = nullptr, *violaver = nullptr, *oobs = nullptr, *orient = nullptr, *odrms = nullptr;
int * index = nullptr, *pair = nullptr, norsel = 0, *orsel = nullptr, *or_label = nullptr;
int nbounds = 0, npairs;
t_inputrec irInstance;
t_inputrec* ir = &irInstance;
init_enxframe(&fr);
- gmx::TopologyInformation topInfo;
+ gmx::TopologyInformation topInfo;
+ std::unique_ptr<gmx_localtop_t> top;
if (!bDisRe)
{
if (bORIRE || bOTEN)
{
{
topInfo.fillFromInputFile(ftp2fn(efTPR, NFILE, fnm));
- gmx_mtop_generate_local_top(*topInfo.mtop(), &top, ir->efep != efepNO);
+ top = std::make_unique<gmx_localtop_t>(topInfo.mtop()->ffparams);
+ gmx_mtop_generate_local_top(*topInfo.mtop(), top.get(), ir->efep != efepNO);
}
- nbounds = get_bounds(&bounds, &index, &pair, &npairs, &top);
+ nbounds = get_bounds(&bounds, &index, &pair, &npairs, top.get());
snew(violaver, npairs);
out_disre = xvgropen(opt2fn("-o", NFILE, fnm), "Sum of Violations", "Time (ps)", "nm", oenv);
xvgr_legend(out_disre, 2, drleg, oenv);
blk_disre = find_block_id_enxframe(&fr, enxDISRE, nullptr);
if (bDisRe && bDRAll && !leg && blk_disre)
{
- t_iatom* fa;
- t_iparams* ip;
-
- fa = top.idef.il[F_DISRES].iatoms;
- ip = top.idef.iparams;
+ const InteractionList& ilist = top->idef.il[F_DISRES];
+ gmx::ArrayRef<const int> fa = ilist.iatoms;
+ const t_iparams* ip = top->idef.iparams.data();
if (blk_disre->nsub != 2 || (blk_disre->sub[0].nr != blk_disre->sub[1].nr))
{
gmx_incons("Number of disre sub-blocks not equal to 2");
}
ndisre = blk_disre->sub[0].nr;
- if (ndisre != top.idef.il[F_DISRES].nr / 3)
+ if (ndisre != ilist.size() / 3)
{
gmx_fatal(FARGS,
"Number of disre pairs in the energy file (%d) does not match the "
"number in the run input file (%d)\n",
- ndisre, top.idef.il[F_DISRES].nr / 3);
+ ndisre, ilist.size() / 3);
}
snew(pairleg, ndisre);
int molb = 0;
gmx::ArrayRef<const t_iparams> iparams,
real massFactorThreshold)
{
- if (molt.ilist[F_CONSTR].size() == 0 && molt.ilist[F_CONSTRNC].size() == 0)
+ if (molt.ilist[F_CONSTR].empty() && molt.ilist[F_CONSTRNC].empty())
{
return false;
}
*/
int ftype_connbond = 0;
int ind_connbond = 0;
- if (molt->ilist[F_CONNBONDS].size() != 0)
+ if (!molt->ilist[F_CONNBONDS].empty())
{
GMX_LOG(logger.info)
.asParagraph()
iloop = gmx_mtop_ilistloop_init(mtop);
while (const InteractionLists* il = gmx_mtop_ilistloop_next(iloop, &nmol))
{
- if (nmol > 1 && (*il)[F_DISRES].size() > 0 && ir->eDisre != edrEnsemble)
+ if (nmol > 1 && !(*il)[F_DISRES].empty() && ir->eDisre != edrEnsemble)
{
gmx_fatal(FARGS,
"NMR distance restraints with multiple copies of the same molecule are "
struct gmx_ffparams_t;
struct gmx_mtop_t;
struct t_forcerec;
-struct t_idef;
struct t_inputrec;
struct gmx_wallcycle;
* stage. Copies the bonded interactions assigned to the GPU
* to device data structures, and updates device buffers that
* may have been updated after search. */
- void updateInteractionListsAndDeviceBuffers(ArrayRef<const int> nbnxnAtomOrder,
- const t_idef& idef,
- void* xqDevice,
- DeviceBuffer<RVec> forceDevice,
- DeviceBuffer<RVec> fshiftDevice);
+ void updateInteractionListsAndDeviceBuffers(ArrayRef<const int> nbnxnAtomOrder,
+ const InteractionDefinitions& idef,
+ void* xqDevice,
+ DeviceBuffer<RVec> forceDevice,
+ DeviceBuffer<RVec> fshiftDevice);
+
/*! \brief Returns whether there are bonded interactions
* assigned to the GPU */
bool haveInteractions() const;
GpuBonded::~GpuBonded() = default;
void GpuBonded::updateInteractionListsAndDeviceBuffers(ArrayRef<const int> /* nbnxnAtomOrder */,
- const t_idef& /* idef */,
+ const InteractionDefinitions& /* idef */,
void* /* xqDevice */,
DeviceBuffer<RVec> /* forceDevice */,
DeviceBuffer<RVec> /* fshiftDevice */)
allocateDeviceBuffer(&d_vTot_, F_NRE, nullptr);
clearDeviceBufferAsync(&d_vTot_, 0, F_NRE, stream_);
- for (int fType = 0; fType < F_NRE; fType++)
- {
- d_iLists_[fType].nr = 0;
- d_iLists_[fType].iatoms = nullptr;
- d_iLists_[fType].nalloc = 0;
- }
-
kernelParams_.d_forceParams = d_forceParams_;
kernelParams_.d_xq = d_xq_;
kernelParams_.d_f = d_f_;
}
//! Return whether function type \p fType in \p idef has perturbed interactions
-static bool fTypeHasPerturbedEntries(const t_idef& idef, int fType)
+static bool fTypeHasPerturbedEntries(const InteractionDefinitions& idef, int fType)
{
GMX_ASSERT(idef.ilsort == ilsortNO_FE || idef.ilsort == ilsortFE_SORTED,
"Perturbed interations should be sorted here");
- const t_ilist& ilist = idef.il[fType];
+ const InteractionList& ilist = idef.il[fType];
- return (idef.ilsort != ilsortNO_FE && idef.numNonperturbedInteractions[fType] != ilist.nr);
+ return (idef.ilsort != ilsortNO_FE && idef.numNonperturbedInteractions[fType] != ilist.size());
}
//! Converts \p src with atom indices in state order to \p dest in nbnxn order
-static void convertIlistToNbnxnOrder(const t_ilist& src,
- HostInteractionList* dest,
- int numAtomsPerInteraction,
- ArrayRef<const int> nbnxnAtomOrder)
+static void convertIlistToNbnxnOrder(const InteractionList& src,
+ HostInteractionList* dest,
+ int numAtomsPerInteraction,
+ ArrayRef<const int> nbnxnAtomOrder)
{
GMX_ASSERT(src.size() == 0 || !nbnxnAtomOrder.empty(), "We need the nbnxn atom order");
* \todo Use DeviceBuffer for the d_xqPtr.
*/
void GpuBonded::Impl::updateInteractionListsAndDeviceBuffers(ArrayRef<const int> nbnxnAtomOrder,
- const t_idef& idef,
- void* d_xqPtr,
- DeviceBuffer<RVec> d_fPtr,
- DeviceBuffer<RVec> d_fShiftPtr)
+ const InteractionDefinitions& idef,
+ void* d_xqPtr,
+ DeviceBuffer<RVec> d_fPtr,
+ DeviceBuffer<RVec> d_fShiftPtr)
{
// TODO wallcycle sub start
haveInteractions_ = false;
* But instead of doing all interactions on the CPU, we can
* still easily handle the types that have no perturbed
* interactions on the GPU. */
- if (idef.il[fType].nr > 0 && !fTypeHasPerturbedEntries(idef, fType))
+ if (!idef.il[fType].empty() && !fTypeHasPerturbedEntries(idef, fType))
{
haveInteractions_ = true;
GpuBonded::~GpuBonded() = default;
-void GpuBonded::updateInteractionListsAndDeviceBuffers(ArrayRef<const int> nbnxnAtomOrder,
- const t_idef& idef,
- void* d_xq,
- DeviceBuffer<RVec> d_f,
- DeviceBuffer<RVec> d_fShift)
+void GpuBonded::updateInteractionListsAndDeviceBuffers(ArrayRef<const int> nbnxnAtomOrder,
+ const InteractionDefinitions& idef,
+ void* d_xq,
+ DeviceBuffer<RVec> d_f,
+ DeviceBuffer<RVec> d_fShift)
{
impl_->updateInteractionListsAndDeviceBuffers(nbnxnAtomOrder, idef, d_xq, d_f, d_fShift);
}
* stage. Copies the bonded interactions assigned to the GPU
* to device data structures, and updates device buffers that
* may have been updated after search. */
- void updateInteractionListsAndDeviceBuffers(ArrayRef<const int> nbnxnAtomOrder,
- const t_idef& idef,
- void* xqDevice,
- DeviceBuffer<RVec> forceDevice,
- DeviceBuffer<RVec> fshiftDevice);
+ void updateInteractionListsAndDeviceBuffers(ArrayRef<const int> nbnxnAtomOrder,
+ const InteractionDefinitions& idef,
+ void* xqDevice,
+ DeviceBuffer<RVec> forceDevice,
+ DeviceBuffer<RVec> fshiftDevice);
/*! \brief Launches bonded kernel on a GPU */
template<bool calcVir, bool calcEner>
//! Tells whether there are any interaction in iLists.
bool haveInteractions_;
//! Interaction lists on the device.
- t_ilist d_iLists_[F_NRE];
+ t_ilist d_iLists_[F_NRE] = {};
//! Bonded parameters for device-side use.
t_iparams* d_forceParams_ = nullptr;
//! Position-charge vector on the device.
/*! \brief Calculate one element of the list of bonded interactions
for this thread */
-real calc_one_bond(int thread,
- int ftype,
- const t_idef* idef,
- ArrayRef<const int> iatoms,
- const int numNonperturbedInteractions,
- const WorkDivision& workDivision,
- const rvec x[],
- rvec4 f[],
- rvec fshift[],
- const t_forcerec* fr,
- const t_pbc* pbc,
- const t_graph* g,
- gmx_grppairener_t* grpp,
- t_nrnb* nrnb,
- const real* lambda,
- real* dvdl,
- const t_mdatoms* md,
- t_fcdata* fcd,
- const gmx::StepWorkload& stepWork,
- int* global_atom_index)
+real calc_one_bond(int thread,
+ int ftype,
+ const InteractionDefinitions& idef,
+ ArrayRef<const int> iatoms,
+ const int numNonperturbedInteractions,
+ const WorkDivision& workDivision,
+ const rvec x[],
+ rvec4 f[],
+ rvec fshift[],
+ const t_forcerec* fr,
+ const t_pbc* pbc,
+ const t_graph* g,
+ gmx_grppairener_t* grpp,
+ t_nrnb* nrnb,
+ const real* lambda,
+ real* dvdl,
+ const t_mdatoms* md,
+ t_fcdata* fcd,
+ const gmx::StepWorkload& stepWork,
+ int* global_atom_index)
{
- GMX_ASSERT(idef->ilsort == ilsortNO_FE || idef->ilsort == ilsortFE_SORTED,
+ GMX_ASSERT(idef.ilsort == ilsortNO_FE || idef.ilsort == ilsortFE_SORTED,
"The topology should be marked either as no FE or sorted on FE");
const bool havePerturbedInteractions =
- (idef->ilsort == ilsortFE_SORTED && numNonperturbedInteractions < iatoms.ssize());
+ (idef.ilsort == ilsortFE_SORTED && numNonperturbedInteractions < iatoms.ssize());
BondedKernelFlavor flavor =
selectBondedKernelFlavor(stepWork, fr->use_simd_kernels, havePerturbedInteractions);
int efptFTYPE;
const int nb0 = workDivision.bound(ftype, thread);
const int nbn = workDivision.bound(ftype, thread + 1) - nb0;
+ ArrayRef<const t_iparams> iparams = idef.iparams;
+
real v = 0;
if (!isPairInteraction(ftype))
{
nice to account to its own subtimer, but first
wallcycle needs to be extended to support calling from
multiple threads. */
- v = cmap_dihs(nbn, iatoms.data() + nb0, idef->iparams, idef->cmap_grid, x, f, fshift,
+ v = cmap_dihs(nbn, iatoms.data() + nb0, iparams.data(), &idef.cmap_grid, x, f, fshift,
pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]), md, fcd, global_atom_index);
}
else
{
- v = calculateSimpleBond(ftype, nbn, iatoms.data() + nb0, idef->iparams, x, f, fshift,
+ v = calculateSimpleBond(ftype, nbn, iatoms.data() + nb0, iparams.data(), x, f, fshift,
pbc, g, lambda[efptFTYPE], &(dvdl[efptFTYPE]), md, fcd,
global_atom_index, flavor);
}
/* TODO The execution time for pairs might be nice to account
to its own subtimer, but first wallcycle needs to be
extended to support calling from multiple threads. */
- do_pairs(ftype, nbn, iatoms.data() + nb0, idef->iparams, x, f, fshift, pbc, g, lambda, dvdl,
- md, fr, havePerturbedInteractions, stepWork, grpp, global_atom_index);
+ do_pairs(ftype, nbn, iatoms.data() + nb0, iparams.data(), x, f, fshift, pbc, g, lambda,
+ dvdl, md, fr, havePerturbedInteractions, stepWork, grpp, global_atom_index);
}
if (thread == 0)
/*! \brief Compute the bonded part of the listed forces, parallelized over threads
*/
-static void calcBondedForces(const t_idef* idef,
- const rvec x[],
- const t_forcerec* fr,
- const t_pbc* pbc_null,
- const t_graph* g,
- rvec* fshiftMasterBuffer,
- gmx_enerdata_t* enerd,
- t_nrnb* nrnb,
- const real* lambda,
- real* dvdl,
- const t_mdatoms* md,
- t_fcdata* fcd,
- const gmx::StepWorkload& stepWork,
- int* global_atom_index)
+static void calcBondedForces(const InteractionDefinitions& idef,
+ const rvec x[],
+ const t_forcerec* fr,
+ const t_pbc* pbc_null,
+ const t_graph* g,
+ rvec* fshiftMasterBuffer,
+ gmx_enerdata_t* enerd,
+ t_nrnb* nrnb,
+ const real* lambda,
+ real* dvdl,
+ const t_mdatoms* md,
+ t_fcdata* fcd,
+ const gmx::StepWorkload& stepWork,
+ int* global_atom_index)
{
bonded_threading_t* bt = fr->bondedThreading;
/* Loop over all bonded force types to calculate the bonded forces */
for (ftype = 0; (ftype < F_NRE); ftype++)
{
- const t_ilist& ilist = idef->il[ftype];
- if (ilist.nr > 0 && ftype_is_bonded_potential(ftype))
+ const InteractionList& ilist = idef.il[ftype];
+ if (!ilist.empty() && ftype_is_bonded_potential(ftype))
{
- ArrayRef<const int> iatoms = gmx::constArrayRefFromArray(ilist.iatoms, ilist.nr);
+ ArrayRef<const int> iatoms = gmx::makeConstArrayRef(ilist.iatoms);
v = calc_one_bond(
- thread, ftype, idef, iatoms, idef->numNonperturbedInteractions[ftype],
+ thread, ftype, idef, iatoms, idef.numNonperturbedInteractions[ftype],
fr->bondedThreading->workDivision, x, ft, fshift, fr, pbc_null, g, grpp,
nrnb, lambda, dvdlt, md, fcd, stepWork, global_atom_index);
epot[ftype] += v;
}
}
-bool haveRestraints(const t_idef& idef, const t_fcdata& fcd)
+bool haveRestraints(const InteractionDefinitions& idef, const t_fcdata& fcd)
{
- return ((idef.il[F_POSRES].nr > 0) || (idef.il[F_FBPOSRES].nr > 0) || fcd.orires.nr > 0
+ return (!idef.il[F_POSRES].empty() || !idef.il[F_FBPOSRES].empty() || fcd.orires.nr > 0
|| fcd.disres.nres > 0);
}
return fr.bondedThreading->haveBondeds;
}
-bool haveCpuListedForces(const t_forcerec& fr, const t_idef& idef, const t_fcdata& fcd)
+bool haveCpuListedForces(const t_forcerec& fr, const InteractionDefinitions& idef, const t_fcdata& fcd)
{
return haveCpuBondeds(fr) || haveRestraints(idef, fcd);
}
-void calc_listed(const t_commrec* cr,
- const gmx_multisim_t* ms,
- struct gmx_wallcycle* wcycle,
- const t_idef* idef,
- const rvec x[],
- history_t* hist,
- gmx::ForceOutputs* forceOutputs,
- const t_forcerec* fr,
- const struct t_pbc* pbc,
- const struct t_pbc* pbc_full,
- const struct t_graph* g,
- gmx_enerdata_t* enerd,
- t_nrnb* nrnb,
- const real* lambda,
- const t_mdatoms* md,
- t_fcdata* fcd,
- int* global_atom_index,
- const gmx::StepWorkload& stepWork)
+void calc_listed(const t_commrec* cr,
+ const gmx_multisim_t* ms,
+ struct gmx_wallcycle* wcycle,
+ const InteractionDefinitions& idef,
+ const rvec x[],
+ history_t* hist,
+ gmx::ForceOutputs* forceOutputs,
+ const t_forcerec* fr,
+ const struct t_pbc* pbc,
+ const struct t_pbc* pbc_full,
+ const struct t_graph* g,
+ gmx_enerdata_t* enerd,
+ t_nrnb* nrnb,
+ const real* lambda,
+ const t_mdatoms* md,
+ t_fcdata* fcd,
+ int* global_atom_index,
+ const gmx::StepWorkload& stepWork)
{
const t_pbc* pbc_null;
bonded_threading_t* bt = fr->bondedThreading;
pbc_null = nullptr;
}
- if (haveRestraints(*idef, *fcd))
+ if (haveRestraints(idef, *fcd))
{
/* TODO Use of restraints triggers further function calls
inside the loop over calc_one_bond(), but those are too
restraints, anyway. */
wallcycle_sub_start(wcycle, ewcsRESTRAINTS);
- if (idef->il[F_POSRES].nr > 0)
+ if (!idef.il[F_POSRES].empty())
{
posres_wrapper(nrnb, idef, pbc_full, x, enerd, lambda, fr, &forceOutputs->forceWithVirial());
}
- if (idef->il[F_FBPOSRES].nr > 0)
+ if (!idef.il[F_FBPOSRES].empty())
{
fbposres_wrapper(nrnb, idef, pbc_full, x, enerd, fr, &forceOutputs->forceWithVirial());
}
*/
GMX_RELEASE_ASSERT(fr->pbcType == PbcType::No || g != nullptr,
"With orientation restraints molecules should be whole");
- enerd->term[F_ORIRESDEV] = calc_orires_dev(ms, idef->il[F_ORIRES].nr, idef->il[F_ORIRES].iatoms,
- idef->iparams, md, x, pbc_null, fcd, hist);
+ enerd->term[F_ORIRESDEV] =
+ calc_orires_dev(ms, idef.il[F_ORIRES].size(), idef.il[F_ORIRES].iatoms.data(),
+ idef.iparams.data(), md, x, pbc_null, fcd, hist);
}
if (fcd->disres.nres > 0)
{
- calc_disres_R_6(cr, ms, idef->il[F_DISRES].nr, idef->il[F_DISRES].iatoms, x, pbc_null,
- fcd, hist);
+ calc_disres_R_6(cr, ms, idef.il[F_DISRES].size(), idef.il[F_DISRES].iatoms.data(), x,
+ pbc_null, fcd, hist);
}
wallcycle_sub_stop(wcycle, ewcsRESTRAINTS);
}
}
-void calc_listed_lambda(const t_idef* idef,
- const rvec x[],
- const t_forcerec* fr,
- const struct t_pbc* pbc,
- const struct t_graph* g,
- gmx_grppairener_t* grpp,
- real* epot,
- gmx::ArrayRef<real> dvdl,
- t_nrnb* nrnb,
- const real* lambda,
- const t_mdatoms* md,
- t_fcdata* fcd,
- int* global_atom_index)
+void calc_listed_lambda(const InteractionDefinitions& idef,
+ const rvec x[],
+ const t_forcerec* fr,
+ const struct t_pbc* pbc,
+ const struct t_graph* g,
+ gmx_grppairener_t* grpp,
+ real* epot,
+ gmx::ArrayRef<real> dvdl,
+ t_nrnb* nrnb,
+ const real* lambda,
+ const t_mdatoms* md,
+ t_fcdata* fcd,
+ int* global_atom_index)
{
real v;
rvec4* f;
rvec* fshift;
const t_pbc* pbc_null;
- t_idef idef_fe;
WorkDivision& workDivision = fr->bondedThreading->foreignLambdaWorkDivision;
if (fr->bMolPBC)
pbc_null = nullptr;
}
- /* Copy the whole idef, so we can modify the contents locally */
- idef_fe = *idef;
-
/* We already have the forces, so we use temp buffers here */
// TODO: Get rid of these allocations by using permanent force buffers
snew(f, fr->natoms_force);
{
if (ftype_is_bonded_potential(ftype))
{
- const t_ilist& ilist = idef->il[ftype];
+ const InteractionList& ilist = idef.il[ftype];
/* Create a temporary iatom list with only perturbed interactions */
- const int numNonperturbed = idef->numNonperturbedInteractions[ftype];
- ArrayRef<const int> iatoms = gmx::constArrayRefFromArray(ilist.iatoms + numNonperturbed,
- ilist.nr - numNonperturbed);
- t_ilist& ilist_fe = idef_fe.il[ftype];
- /* Set the work range of thread 0 to the perturbed bondeds */
- workDivision.setBound(ftype, 0, 0);
- workDivision.setBound(ftype, 1, iatoms.ssize());
-
- if (ilist_fe.nr > 0)
+ const int numNonperturbed = idef.numNonperturbedInteractions[ftype];
+ ArrayRef<const int> iatomsPerturbed = gmx::constArrayRefFromArray(
+ ilist.iatoms.data() + numNonperturbed, ilist.size() - numNonperturbed);
+ if (!iatomsPerturbed.empty())
{
+ /* Set the work range of thread 0 to the perturbed bondeds */
+ workDivision.setBound(ftype, 0, 0);
+ workDivision.setBound(ftype, 1, iatomsPerturbed.ssize());
+
gmx::StepWorkload tempFlags;
tempFlags.computeEnergy = true;
- v = calc_one_bond(0, ftype, idef, iatoms, iatoms.ssize(), workDivision, x, f,
- fshift, fr, pbc_null, g, grpp, nrnb, lambda, dvdl.data(), md, fcd,
- tempFlags, global_atom_index);
+ v = calc_one_bond(0, ftype, idef, iatomsPerturbed, iatomsPerturbed.ssize(),
+ workDivision, x, f, fshift, fr, pbc_null, g, grpp, nrnb, lambda,
+ dvdl.data(), md, fcd, tempFlags, global_atom_index);
epot[ftype] += v;
}
}
sfree(f);
}
-void do_force_listed(struct gmx_wallcycle* wcycle,
- const matrix box,
- const t_lambda* fepvals,
- const t_commrec* cr,
- const gmx_multisim_t* ms,
- const t_idef* idef,
- const rvec x[],
- history_t* hist,
- gmx::ForceOutputs* forceOutputs,
- const t_forcerec* fr,
- const struct t_pbc* pbc,
- const struct t_graph* graph,
- gmx_enerdata_t* enerd,
- t_nrnb* nrnb,
- const real* lambda,
- const t_mdatoms* md,
- t_fcdata* fcd,
- int* global_atom_index,
- const gmx::StepWorkload& stepWork)
+void do_force_listed(struct gmx_wallcycle* wcycle,
+ const matrix box,
+ const t_lambda* fepvals,
+ const t_commrec* cr,
+ const gmx_multisim_t* ms,
+ const InteractionDefinitions& idef,
+ const rvec x[],
+ history_t* hist,
+ gmx::ForceOutputs* forceOutputs,
+ const t_forcerec* fr,
+ const struct t_pbc* pbc,
+ const struct t_graph* graph,
+ gmx_enerdata_t* enerd,
+ t_nrnb* nrnb,
+ const real* lambda,
+ const t_mdatoms* md,
+ t_fcdata* fcd,
+ int* global_atom_index,
+ const gmx::StepWorkload& stepWork)
{
t_pbc pbc_full; /* Full PBC is needed for position restraints */
return;
}
- if ((idef->il[F_POSRES].nr > 0) || (idef->il[F_FBPOSRES].nr > 0))
+ if (!idef.il[F_POSRES].empty() || !idef.il[F_FBPOSRES].empty())
{
/* Not enough flops to bother counting */
set_pbc(&pbc_full, fr->pbcType, box);
real dvdl[efptNR] = { 0 };
posres_wrapper_lambda(wcycle, fepvals, idef, &pbc_full, x, enerd, lambda, fr);
- if (idef->ilsort != ilsortNO_FE)
+ if (idef.ilsort != ilsortNO_FE)
{
wallcycle_sub_start(wcycle, ewcsLISTED_FEP);
- if (idef->ilsort != ilsortFE_SORTED)
+ if (idef.ilsort != ilsortFE_SORTED)
{
gmx_incons("The bonded interactions are not sorted for free energy");
}
struct gmx_grppairener_t;
struct gmx_multisim_t;
class history_t;
+class InteractionDefinitions;
struct t_commrec;
struct t_fcdata;
struct t_forcerec;
-struct t_idef;
struct t_graph;
struct t_lambda;
struct t_mdatoms;
*
* Note that pbc_full is used only for position restraints, and is
* not initialized if there are none. */
-void calc_listed(const t_commrec* cr,
- const gmx_multisim_t* ms,
- struct gmx_wallcycle* wcycle,
- const t_idef* idef,
- const rvec x[],
- history_t* hist,
- gmx::ForceOutputs* forceOutputs,
- const t_forcerec* fr,
- const struct t_pbc* pbc,
- const struct t_pbc* pbc_full,
- const struct t_graph* g,
- gmx_enerdata_t* enerd,
- t_nrnb* nrnb,
- const real* lambda,
- const t_mdatoms* md,
- struct t_fcdata* fcd,
- int* ddgatindex,
- const gmx::StepWorkload& stepWork);
+void calc_listed(const t_commrec* cr,
+ const gmx_multisim_t* ms,
+ struct gmx_wallcycle* wcycle,
+ const InteractionDefinitions& idef,
+ const rvec x[],
+ history_t* hist,
+ gmx::ForceOutputs* forceOutputs,
+ const t_forcerec* fr,
+ const struct t_pbc* pbc,
+ const struct t_pbc* pbc_full,
+ const struct t_graph* g,
+ gmx_enerdata_t* enerd,
+ t_nrnb* nrnb,
+ const real* lambda,
+ const t_mdatoms* md,
+ struct t_fcdata* fcd,
+ int* ddgatindex,
+ const gmx::StepWorkload& stepWork);
/*! \brief As calc_listed(), but only determines the potential energy
* for the perturbed interactions.
*
* The shift forces in fr are not affected. */
-void calc_listed_lambda(const t_idef* idef,
- const rvec x[],
- const t_forcerec* fr,
- const struct t_pbc* pbc,
- const struct t_graph* g,
- gmx_grppairener_t* grpp,
- real* epot,
- gmx::ArrayRef<real> dvdl,
- t_nrnb* nrnb,
- const real* lambda,
- const t_mdatoms* md,
- struct t_fcdata* fcd,
- int* global_atom_index);
+void calc_listed_lambda(const InteractionDefinitions& idef,
+ const rvec x[],
+ const t_forcerec* fr,
+ const struct t_pbc* pbc,
+ const struct t_graph* g,
+ gmx_grppairener_t* grpp,
+ real* epot,
+ gmx::ArrayRef<real> dvdl,
+ t_nrnb* nrnb,
+ const real* lambda,
+ const t_mdatoms* md,
+ struct t_fcdata* fcd,
+ int* global_atom_index);
/*! \brief Do all aspects of energy and force calculations for mdrun
* on the set of listed interactions */
-void do_force_listed(struct gmx_wallcycle* wcycle,
- const matrix box,
- const t_lambda* fepvals,
- const t_commrec* cr,
- const gmx_multisim_t* ms,
- const t_idef* idef,
- const rvec x[],
- history_t* hist,
- gmx::ForceOutputs* forceOutputs,
- const t_forcerec* fr,
- const struct t_pbc* pbc,
- const struct t_graph* graph,
- gmx_enerdata_t* enerd,
- t_nrnb* nrnb,
- const real* lambda,
- const t_mdatoms* md,
- struct t_fcdata* fcd,
- int* global_atom_index,
- const gmx::StepWorkload& stepWork);
+void do_force_listed(struct gmx_wallcycle* wcycle,
+ const matrix box,
+ const t_lambda* fepvals,
+ const t_commrec* cr,
+ const gmx_multisim_t* ms,
+ const InteractionDefinitions& idef,
+ const rvec x[],
+ history_t* hist,
+ gmx::ForceOutputs* forceOutputs,
+ const t_forcerec* fr,
+ const struct t_pbc* pbc,
+ const struct t_graph* graph,
+ gmx_enerdata_t* enerd,
+ t_nrnb* nrnb,
+ const real* lambda,
+ const t_mdatoms* md,
+ struct t_fcdata* fcd,
+ int* global_atom_index,
+ const gmx::StepWorkload& stepWork);
/*! \brief Returns true if there are position, distance or orientation restraints. */
-bool haveRestraints(const t_idef& idef, const t_fcdata& fcd);
+bool haveRestraints(const InteractionDefinitions& idef, const t_fcdata& fcd);
/*! \brief Returns true if there are CPU (i.e. not GPU-offloaded) bonded interactions to compute. */
bool haveCpuBondeds(const t_forcerec& fr);
* NOTE: the current implementation returns true if there are position restraints
* or any bonded interactions computed on the CPU.
*/
-bool haveCpuListedForces(const t_forcerec& fr, const t_idef& idef, const t_fcdata& fcd);
+bool haveCpuListedForces(const t_forcerec& fr, const InteractionDefinitions& idef, const t_fcdata& fcd);
#endif
/*! \brief struct for passing all data required for a function type */
typedef struct
{
- const t_ilist* il; /**< pointer to t_ilist entry corresponding to ftype */
- int ftype; /**< the function type index */
- int nat; /**< nr of atoms involved in a single ftype interaction */
+ const InteractionList* il; /**< pointer to t_ilist entry corresponding to ftype */
+ int ftype; /**< the function type index */
+ int nat; /**< nr of atoms involved in a single ftype interaction */
} ilist_data_t;
/*! \brief Divides listed interactions over threads
for (f = 0; f < numType; f++)
{
/* Sum #bondeds*#atoms_per_bond over all bonded types */
- nat_tot += ild[f].il->nr / (ild[f].nat + 1) * ild[f].nat;
+ nat_tot += ild[f].il->size() / (ild[f].nat + 1) * ild[f].nat;
/* The start bound for thread 0 is 0 for all interactions */
ind[f] = 0;
/* Initialize the next atom index array */
- assert(ild[f].il->nr > 0);
+ assert(!ild[f].il->empty());
at_ind[f] = ild[f].il->iatoms[1];
}
nat_sum += ild[f_min].nat;
/* Update the first unassigned atom index for this type */
- if (ind[f_min] < ild[f_min].il->nr)
+ if (ind[f_min] < ild[f_min].il->size())
{
at_ind[f_min] = ild[f_min].il->iatoms[ind[f_min] + 1];
}
for (f = 0; f < numType; f++)
{
- assert(ind[f] == ild[f].il->nr);
+ assert(ind[f] == ild[f].il->size());
}
}
//! Return whether function type \p ftype in \p idef has perturbed interactions
-static bool ftypeHasPerturbedEntries(const t_idef& idef, int ftype)
+static bool ftypeHasPerturbedEntries(const InteractionDefinitions& idef, int ftype)
{
GMX_ASSERT(idef.ilsort == ilsortNO_FE || idef.ilsort == ilsortFE_SORTED,
"Perturbed interations should be sorted here");
- const t_ilist& ilist = idef.il[ftype];
+ const InteractionList& ilist = idef.il[ftype];
- return (idef.ilsort != ilsortNO_FE && idef.numNonperturbedInteractions[ftype] != ilist.nr);
+ return (idef.ilsort != ilsortNO_FE && idef.numNonperturbedInteractions[ftype] != ilist.size());
}
//! Divides bonded interactions over threads and GPU
-static void divide_bondeds_over_threads(bonded_threading_t* bt, bool useGpuForBondeds, const t_idef& idef)
+static void divide_bondeds_over_threads(bonded_threading_t* bt,
+ bool useGpuForBondeds,
+ const InteractionDefinitions& idef)
{
ilist_data_t ild[F_NRE];
GMX_ASSERT(bt->nthreads > 0, "Must have positive number of threads");
const int numThreads = bt->nthreads;
+ gmx::ArrayRef<const t_iparams> iparams = idef.iparams;
+
bt->haveBondeds = false;
int numType = 0;
size_t fTypeGpuIndex = 0;
continue;
}
- const t_ilist& il = idef.il[fType];
- int nrToAssignToCpuThreads = il.nr;
+ const InteractionList& il = idef.il[fType];
+ int nrToAssignToCpuThreads = il.size();
if (useGpuForBondeds && fTypeGpuIndex < gmx::fTypesOnGpu.size()
&& gmx::fTypesOnGpu[fTypeGpuIndex] == fType)
* end up on the same thread.
*/
while (nr_t > 0 && nr_t < nrToAssignToCpuThreads
- && idef.iparams[il.iatoms[nr_t]].disres.label
- == idef.iparams[il.iatoms[nr_t - stride]].disres.label)
+ && iparams[il.iatoms[nr_t]].disres.label
+ == iparams[il.iatoms[nr_t - stride]].disres.label)
{
nr_t += stride;
}
fprintf(debug, "Division of bondeds over threads:\n");
for (f = 0; f < F_NRE; f++)
{
- if (ftype_is_bonded_potential(f) && idef.il[f].nr > 0)
+ if (ftype_is_bonded_potential(f) && !idef.il[f].empty())
{
int t;
}
//! Construct a reduction mask for which parts (blocks) of the force array are touched on which thread task
-static void calc_bonded_reduction_mask(int natoms,
- f_thread_t* f_thread,
- const t_idef& idef,
- int thread,
- const bonded_threading_t& bondedThreading)
+static void calc_bonded_reduction_mask(int natoms,
+ f_thread_t* f_thread,
+ const InteractionDefinitions& idef,
+ int thread,
+ const bonded_threading_t& bondedThreading)
{
static_assert(BITMASK_SIZE == GMX_OPENMP_MAX_THREADS,
"For the error message below we assume these two are equal.");
{
if (ftype_is_bonded_potential(ftype))
{
- int nb = idef.il[ftype].nr;
+ int nb = idef.il[ftype].size();
if (nb > 0)
{
int nat1 = interaction_function[ftype].nratoms + 1;
}
}
-void setup_bonded_threading(bonded_threading_t* bt, int numAtoms, bool useGpuForBondeds, const t_idef& idef)
+void setup_bonded_threading(bonded_threading_t* bt,
+ int numAtoms,
+ bool useGpuForBondeds,
+ const InteractionDefinitions& idef)
{
int ctot = 0;
#include <cstdio>
struct bonded_threading_t;
-struct t_idef;
+class InteractionDefinitions;
/*! \brief Divide the listed interactions over the threads and GPU
*
* This should be called each time the bonded setup changes;
* i.e. at start-up without domain decomposition and at DD.
*/
-void setup_bonded_threading(bonded_threading_t* bt, int numAtoms, bool useGpuForBondes, const t_idef& idef);
+void setup_bonded_threading(bonded_threading_t* bt,
+ int numAtoms,
+ bool useGpuForBondes,
+ const InteractionDefinitions& idef);
//! Destructor.
void tear_down_bonded_threading(bonded_threading_t* bt);
int nmol;
while (const InteractionLists* il = gmx_mtop_ilistloop_next(iloop, &nmol))
{
- if (nmol > 1 && (*il)[F_ORIRES].size() > 0)
+ const int numOrires = (*il)[F_ORIRES].size();
+ if (nmol > 1 && numOrires > 0)
{
gmx_fatal(FARGS,
"Found %d copies of a molecule with orientation restrains while the current "
nmol);
}
- for (int i = 0; i < (*il)[F_ORIRES].size(); i += 3)
+ for (int i = 0; i < numOrires; i += 3)
{
int type = (*il)[F_ORIRES].iatoms[i];
int ex = mtop->ffparams.iparams[type].orires.ex;
} // namespace
-void posres_wrapper(t_nrnb* nrnb,
- const t_idef* idef,
- const struct t_pbc* pbc,
- const rvec* x,
- gmx_enerdata_t* enerd,
- const real* lambda,
- const t_forcerec* fr,
- gmx::ForceWithVirial* forceWithVirial)
+void posres_wrapper(t_nrnb* nrnb,
+ const InteractionDefinitions& idef,
+ const struct t_pbc* pbc,
+ const rvec* x,
+ gmx_enerdata_t* enerd,
+ const real* lambda,
+ const t_forcerec* fr,
+ gmx::ForceWithVirial* forceWithVirial)
{
real v, dvdl;
dvdl = 0;
- v = posres<true>(idef->il[F_POSRES].nr, idef->il[F_POSRES].iatoms, idef->iparams_posres, x,
- forceWithVirial, fr->pbcType == PbcType::No ? nullptr : pbc, lambda[efptRESTRAINT],
- &dvdl, fr->rc_scaling, fr->pbcType, fr->posres_com, fr->posres_comB);
+ v = posres<true>(idef.il[F_POSRES].size(), idef.il[F_POSRES].iatoms.data(),
+ idef.iparams_posres.data(), x, forceWithVirial,
+ fr->pbcType == PbcType::No ? nullptr : pbc, lambda[efptRESTRAINT], &dvdl,
+ fr->rc_scaling, fr->pbcType, fr->posres_com, fr->posres_comB);
enerd->term[F_POSRES] += v;
/* If just the force constant changes, the FEP term is linear,
* but if k changes, it is not.
*/
enerd->dvdl_nonlin[efptRESTRAINT] += dvdl;
- inc_nrnb(nrnb, eNR_POSRES, gmx::exactDiv(idef->il[F_POSRES].nr, 2));
+ inc_nrnb(nrnb, eNR_POSRES, gmx::exactDiv(idef.il[F_POSRES].size(), 2));
}
-void posres_wrapper_lambda(struct gmx_wallcycle* wcycle,
- const t_lambda* fepvals,
- const t_idef* idef,
- const struct t_pbc* pbc,
- const rvec x[],
- gmx_enerdata_t* enerd,
- const real* lambda,
- const t_forcerec* fr)
+void posres_wrapper_lambda(struct gmx_wallcycle* wcycle,
+ const t_lambda* fepvals,
+ const InteractionDefinitions& idef,
+ const struct t_pbc* pbc,
+ const rvec x[],
+ gmx_enerdata_t* enerd,
+ const real* lambda,
+ const t_forcerec* fr)
{
real v;
- if (0 == idef->il[F_POSRES].nr)
+ if (idef.il[F_POSRES].empty())
{
return;
}
const real lambda_dum =
(i == 0 ? lambda[efptRESTRAINT] : fepvals->all_lambda[efptRESTRAINT][i - 1]);
- v = posres<false>(idef->il[F_POSRES].nr, idef->il[F_POSRES].iatoms, idef->iparams_posres, x,
- nullptr, fr->pbcType == PbcType::No ? nullptr : pbc, lambda_dum, &dvdl,
+ v = posres<false>(idef.il[F_POSRES].size(), idef.il[F_POSRES].iatoms.data(),
+ idef.iparams_posres.data(), x, nullptr,
+ fr->pbcType == PbcType::No ? nullptr : pbc, lambda_dum, &dvdl,
fr->rc_scaling, fr->pbcType, fr->posres_com, fr->posres_comB);
enerd->enerpart_lambda[i] += v;
enerd->dhdlLambda[i] += dvdl;
/*! \brief Helper function that wraps calls to fbposres for
free-energy perturbation */
-void fbposres_wrapper(t_nrnb* nrnb,
- const t_idef* idef,
- const struct t_pbc* pbc,
- const rvec* x,
- gmx_enerdata_t* enerd,
- const t_forcerec* fr,
- gmx::ForceWithVirial* forceWithVirial)
+void fbposres_wrapper(t_nrnb* nrnb,
+ const InteractionDefinitions& idef,
+ const struct t_pbc* pbc,
+ const rvec* x,
+ gmx_enerdata_t* enerd,
+ const t_forcerec* fr,
+ gmx::ForceWithVirial* forceWithVirial)
{
real v;
- v = fbposres(idef->il[F_FBPOSRES].nr, idef->il[F_FBPOSRES].iatoms, idef->iparams_fbposres, x,
- forceWithVirial, fr->pbcType == PbcType::No ? nullptr : pbc, fr->rc_scaling,
- fr->pbcType, fr->posres_com);
+ v = fbposres(idef.il[F_FBPOSRES].size(), idef.il[F_FBPOSRES].iatoms.data(),
+ idef.iparams_fbposres.data(), x, forceWithVirial,
+ fr->pbcType == PbcType::No ? nullptr : pbc, fr->rc_scaling, fr->pbcType, fr->posres_com);
enerd->term[F_FBPOSRES] += v;
- inc_nrnb(nrnb, eNR_FBPOSRES, gmx::exactDiv(idef->il[F_FBPOSRES].nr, 2));
+ inc_nrnb(nrnb, eNR_FBPOSRES, gmx::exactDiv(idef.il[F_FBPOSRES].size(), 2));
}
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2013,2014,2015,2017,2019, by the GROMACS development team, led by
+ * Copyright (c) 2013,2014,2015,2017,2019,2020, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
struct gmx_enerdata_t;
struct gmx_wallcycle;
struct t_forcerec;
-struct t_idef;
+class InteractionDefinitions;
struct t_lambda;
struct t_nrnb;
struct t_pbc;
}
/*! \brief Helper function that wraps calls to posres */
-void posres_wrapper(t_nrnb* nrnb,
- const t_idef* idef,
- const struct t_pbc* pbc,
- const rvec* x,
- gmx_enerdata_t* enerd,
- const real* lambda,
- const t_forcerec* fr,
- gmx::ForceWithVirial* forceWithVirial);
+void posres_wrapper(t_nrnb* nrnb,
+ const InteractionDefinitions& idef,
+ const struct t_pbc* pbc,
+ const rvec* x,
+ gmx_enerdata_t* enerd,
+ const real* lambda,
+ const t_forcerec* fr,
+ gmx::ForceWithVirial* forceWithVirial);
/*! \brief Helper function that wraps calls to posres for free-energy
pertubation */
-void posres_wrapper_lambda(struct gmx_wallcycle* wcycle,
- const t_lambda* fepvals,
- const t_idef* idef,
- const struct t_pbc* pbc,
- const rvec x[],
- gmx_enerdata_t* enerd,
- const real* lambda,
- const t_forcerec* fr);
+void posres_wrapper_lambda(struct gmx_wallcycle* wcycle,
+ const t_lambda* fepvals,
+ const InteractionDefinitions& idef,
+ const struct t_pbc* pbc,
+ const rvec x[],
+ gmx_enerdata_t* enerd,
+ const real* lambda,
+ const t_forcerec* fr);
/*! \brief Helper function that wraps calls to fbposres for
free-energy perturbation */
-void fbposres_wrapper(t_nrnb* nrnb,
- const t_idef* idef,
- const struct t_pbc* pbc,
- const rvec* x,
- gmx_enerdata_t* enerd,
- const t_forcerec* fr,
- gmx::ForceWithVirial* forceWithVirial);
+void fbposres_wrapper(t_nrnb* nrnb,
+ const InteractionDefinitions& idef,
+ const struct t_pbc* pbc,
+ const rvec* x,
+ gmx_enerdata_t* enerd,
+ const t_forcerec* fr,
+ gmx::ForceWithVirial* forceWithVirial);
#endif
* This file is part of the GROMACS molecular simulation package.
*
* Copyright (c) 2012-2018, The GROMACS development team.
- * Copyright (c) 2019, by the GROMACS development team, led by
+ * Copyright (c) 2019,2020, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
for (const int atom : block)
{
const auto& ilist = moltype.ilist[F_SETTLE];
- GMX_RELEASE_ASSERT(ilist.size() > 0,
+ GMX_RELEASE_ASSERT(!ilist.empty(),
"There should be at least one settle in this moltype");
for (int i = 0; i < ilist.size(); i += 1 + NRAL(F_SETTLE))
{
int numConstraints,
int numSettles);
~Impl();
- void setConstraints(const gmx_localtop_t& top, const t_mdatoms& md);
+ void setConstraints(gmx_localtop_t* top, const t_mdatoms& md);
bool apply(bool bLog,
bool bEner,
int64_t step,
//! Pointer to the global topology - only used for printing warnings.
const gmx_mtop_t& mtop;
//! Parameters for the interactions in this domain.
- const t_idef* idef = nullptr;
+ const InteractionDefinitions* idef = nullptr;
//! Data about atoms in this domain.
const t_mdatoms& md;
//! Whether we need to do pbc for handling bonds.
{
clear_mat(vir_r_m_dr);
}
- const t_ilist* settle = &idef->il[F_SETTLE];
- nsettle = settle->nr / (1 + NRAL(F_SETTLE));
+ const InteractionList& settle = idef->il[F_SETTLE];
+ nsettle = settle.size() / (1 + NRAL(F_SETTLE));
if (nsettle > 0)
{
if (start_th >= 0 && end_th - start_th > 0)
{
settle_proj(settled, econq, end_th - start_th,
- settle->iatoms + start_th * (1 + NRAL(F_SETTLE)), pbc_null,
+ settle.iatoms.data() + start_th * (1 + NRAL(F_SETTLE)), pbc_null,
x.unpaddedArrayRef(), xprime.unpaddedArrayRef(), min_proj,
calcvir_atom_end, th == 0 ? vir_r_m_dr : vir_r_m_dr_th[th]);
}
*
* \returns a block struct with all constraints for each atom
*/
-template<typename T>
-static ListOfLists<int> makeAtomsToConstraintsList(int numAtoms,
- const T* ilists,
- const t_iparams* iparams,
+static ListOfLists<int> makeAtomsToConstraintsList(int numAtoms,
+ ArrayRef<const InteractionList> ilists,
+ ArrayRef<const t_iparams> iparams,
FlexibleConstraintTreatment flexibleConstraintTreatment)
{
- GMX_ASSERT(flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || iparams != nullptr,
+ GMX_ASSERT(flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || !iparams.empty(),
"With flexible constraint detection we need valid iparams");
std::vector<int> count(numAtoms);
for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
{
- const T& ilist = ilists[ftype];
- const int stride = 1 + NRAL(ftype);
+ const InteractionList& ilist = ilists[ftype];
+ const int stride = 1 + NRAL(ftype);
for (int i = 0; i < ilist.size(); i += stride)
{
if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include
int numConstraints = 0;
for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
{
- const T& ilist = ilists[ftype];
- const int stride = 1 + NRAL(ftype);
+ const InteractionList& ilist = ilists[ftype];
+ const int stride = 1 + NRAL(ftype);
for (int i = 0; i < ilist.size(); i += stride)
{
if (flexibleConstraintTreatment == FlexibleConstraintTreatment::Include
return ListOfLists<int>(std::move(listRanges), std::move(elements));
}
-ListOfLists<int> make_at2con(int numAtoms,
- const t_ilist* ilist,
- const t_iparams* iparams,
- FlexibleConstraintTreatment flexibleConstraintTreatment)
+ListOfLists<int> make_at2con(int numAtoms,
+ ArrayRef<const InteractionList> ilist,
+ ArrayRef<const t_iparams> iparams,
+ FlexibleConstraintTreatment flexibleConstraintTreatment)
{
return makeAtomsToConstraintsList(numAtoms, ilist, iparams, flexibleConstraintTreatment);
}
gmx::ArrayRef<const t_iparams> iparams,
FlexibleConstraintTreatment flexibleConstraintTreatment)
{
- return makeAtomsToConstraintsList(moltype.atoms.nr, moltype.ilist.data(), iparams.data(),
+ return makeAtomsToConstraintsList(moltype.atoms.nr, makeConstArrayRef(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 countFlexibleConstraints(ArrayRef<const InteractionList> ilist, ArrayRef<const t_iparams> iparams)
{
int nflexcon = 0;
for (int ftype = F_CONSTR; ftype <= F_CONSTRNC; ftype++)
return nflexcon;
}
-int countFlexibleConstraints(const t_ilist* ilist, const t_iparams* iparams)
-{
- return countFlexibleConstraintsTemplate(ilist, iparams);
-}
-
//! Returns the index of the settle to which each atom belongs.
static std::vector<int> make_at2settle(int natoms, const InteractionList& ilist)
{
return at2s;
}
-void Constraints::Impl::setConstraints(const gmx_localtop_t& top, const t_mdatoms& md)
+void Constraints::Impl::setConstraints(gmx_localtop_t* top, const t_mdatoms& md)
{
- idef = &top.idef;
+ idef = &top->idef;
if (ncon_tot > 0)
{
*/
if (ir.eConstrAlg == econtLINCS)
{
- set_lincs(top.idef, md, EI_DYNAMICS(ir.eI), cr, lincsd);
+ set_lincs(*idef, md, EI_DYNAMICS(ir.eI), cr, lincsd);
}
if (ir.eConstrAlg == econtSHAKE)
{
{
// We are using the local topology, so there are only
// F_CONSTR constraints.
- make_shake_sblock_dd(shaked, &idef->il[F_CONSTR], cr->dd);
+ make_shake_sblock_dd(shaked, idef->il[F_CONSTR], cr->dd);
}
else
{
- make_shake_sblock_serial(shaked, idef, md);
+ make_shake_sblock_serial(shaked, &top->idef, md);
}
}
}
if (settled)
{
- settle_set_constraints(settled, &idef->il[F_SETTLE], md);
+ settle_set_constraints(settled, idef->il[F_SETTLE], md);
}
/* Make a selection of the local atoms for essential dynamics */
}
}
-void Constraints::setConstraints(const gmx_localtop_t& top, const t_mdatoms& md)
+void Constraints::setConstraints(gmx_localtop_t* top, const t_mdatoms& md)
{
impl_->setConstraints(top, md);
}
for (const gmx_molblock_t& molblock : mtop.molblock)
{
- int count = countFlexibleConstraintsTemplate(mtop.moltype[molblock.type].ilist.data(),
- mtop.ffparams.iparams.data());
+ int count = countFlexibleConstraints(mtop.moltype[molblock.type].ilist, mtop.ffparams.iparams);
nflexcon += molblock.nmol * count;
}
*
* \todo Make this a callback that is called automatically
* once a new domain has been made. */
- void setConstraints(const gmx_localtop_t& top, const t_mdatoms& md);
+ void setConstraints(gmx_localtop_t* top, const t_mdatoms& md);
/*! \brief Applies constraints to coordinates.
*
[[noreturn]] void too_many_constraint_warnings(int eConstrAlg, int warncount);
/*! \brief Returns whether constraint with parameter \p iparamsIndex is a flexible constraint */
-static inline bool isConstraintFlexible(const t_iparams* iparams, int iparamsIndex)
+static inline bool isConstraintFlexible(ArrayRef<const t_iparams> iparams, int iparamsIndex)
{
- GMX_ASSERT(iparams != nullptr, "Need a valid iparams array");
-
return (iparams[iparamsIndex].constr.dA == 0 && iparams[iparamsIndex].constr.dB == 0);
};
*
* \returns a ListOfLists object with all constraints for each atom
*/
-ListOfLists<int> make_at2con(int numAtoms,
- const t_ilist* ilist,
- const t_iparams* iparams,
- FlexibleConstraintTreatment flexibleConstraintTreatment);
+ListOfLists<int> make_at2con(int numAtoms,
+ ArrayRef<const InteractionList> ilist,
+ ArrayRef<const t_iparams> iparams,
+ FlexibleConstraintTreatment flexibleConstraintTreatment);
//! Return the number of flexible constraints in the \c ilist and \c iparams.
-int countFlexibleConstraints(const t_ilist* ilist, const t_iparams* iparams);
+int countFlexibleConstraints(ArrayRef<const InteractionList> ilist, ArrayRef<const t_iparams> iparams);
/*! \brief Returns the constraint iatoms for a constraint number con
* which comes from a list where F_CONSTR and F_CONSTRNC constraints
real r0, r1, r2maxA, r2maxB, rmax, lam0, lam1;
- if (molt->ilist[F_CONSTR].size() == 0 && molt->ilist[F_CONSTRNC].size() == 0)
+ if (molt->ilist[F_CONSTR].empty() && molt->ilist[F_CONSTRNC].empty())
{
return 0;
}
void do_force_lowlevel(t_forcerec* fr,
const t_inputrec* ir,
- const t_idef* idef,
+ const InteractionDefinitions& idef,
const t_commrec* cr,
const gmx_multisim_t* ms,
t_nrnb* nrnb,
/* Check whether we need to take into account PBC in listed interactions. */
const auto needPbcForListedForces =
- fr->bMolPBC && stepWork.computeListedForces && haveCpuListedForces(*fr, *idef, *fcd);
+ fr->bMolPBC && stepWork.computeListedForces && haveCpuListedForces(*fr, idef, *fcd);
if (needPbcForListedForces)
{
/* Since all atoms are in the rectangular or triclinic unit-cell,
struct gmx_vsite_t;
struct gmx_wallcycle;
class history_t;
+class InteractionDefinitions;
struct pull_t;
struct t_commrec;
struct t_fcdata;
struct t_forcerec;
struct t_graph;
-struct t_idef;
struct t_inputrec;
struct t_lambda;
struct t_mdatoms;
void do_force_lowlevel(t_forcerec* fr,
const t_inputrec* ir,
- const t_idef* idef,
+ const InteractionDefinitions& idef,
const t_commrec* cr,
const gmx_multisim_t* ms,
t_nrnb* nrnb,
/*! \brief Check if constraint with topology index constraint_index is connected
* to other constraints, and if so add those connected constraints to our task. */
-static void check_assign_connected(Lincs* li,
- const t_iatom* iatom,
- const t_idef& idef,
- bool bDynamics,
- int a1,
- int a2,
- const ListOfLists<int>& at2con)
+static void check_assign_connected(Lincs* li,
+ gmx::ArrayRef<const int> iatom,
+ const InteractionDefinitions& idef,
+ bool bDynamics,
+ int a1,
+ int a2,
+ const ListOfLists<int>& at2con)
{
/* Currently this function only supports constraint groups
* in which all constraints share at least one atom
/*! \brief Check if constraint with topology index constraint_index is involved
* in a constraint triangle, and if so add the other two constraints
* in the triangle to our task. */
-static void check_assign_triangle(Lincs* li,
- const t_iatom* iatom,
- const t_idef& idef,
- bool bDynamics,
- int constraint_index,
- int a1,
- int a2,
- const ListOfLists<int>& at2con)
+static void check_assign_triangle(Lincs* li,
+ gmx::ArrayRef<const int> iatom,
+ const InteractionDefinitions& idef,
+ bool bDynamics,
+ int constraint_index,
+ int a1,
+ int a2,
+ const ListOfLists<int>& at2con)
{
int nca, cc[32], ca[32];
int c_triangle[2] = { -1, -1 };
}
}
-void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_commrec* cr, Lincs* li)
+void set_lincs(const InteractionDefinitions& idef, const t_mdatoms& md, bool bDynamics, const t_commrec* cr, Lincs* li)
{
- int natoms;
- t_iatom* iatom;
-
li->nc_real = 0;
li->nc = 0;
li->ncc = 0;
}
/* This is the local topology, so there are only F_CONSTR constraints */
- if (idef.il[F_CONSTR].nr == 0)
+ if (idef.il[F_CONSTR].empty())
{
/* There are no constraints,
* we do not need to fill any data structures.
fprintf(debug, "Building the LINCS connectivity\n");
}
+ int natoms;
if (DOMAINDECOMP(cr))
{
if (cr->dd->constraints)
const ListOfLists<int> at2con =
make_at2con(natoms, idef.il, idef.iparams, flexibleConstraintTreatment(bDynamics));
- const int ncon_tot = idef.il[F_CONSTR].nr / 3;
+ const int ncon_tot = idef.il[F_CONSTR].size() / 3;
/* Ensure we have enough padding for aligned loads for each thread */
const int numEntries = ncon_tot + li->ntask * simd_width;
li->tmp4.resize(numEntries);
li->mlambda.resize(numEntries);
- iatom = idef.il[F_CONSTR].iatoms;
+ gmx::ArrayRef<const int> iatom = idef.il[F_CONSTR].iatoms;
li->blnr[0] = li->ncc;
*/
li_task->b0 = li->nc;
+ gmx::ArrayRef<const t_iparams> iparams = idef.iparams;
+
while (con < ncon_tot && li->nc - li_task->b0 < ncon_target)
{
if (li->con_index[con] == -1)
type = iatom[3 * con];
a1 = iatom[3 * con + 1];
a2 = iatom[3 * con + 2];
- lenA = idef.iparams[type].constr.dA;
- lenB = idef.iparams[type].constr.dB;
+ lenA = iparams[type].constr.dA;
+ lenB = iparams[type].constr.dB;
/* Skip the flexible constraints when not doing dynamics */
if (bDynamics || lenA != 0 || lenB != 0)
{
struct gmx_mtop_t;
struct gmx_multisim_t;
+class InteractionDefinitions;
struct t_commrec;
-struct t_idef;
struct t_inputrec;
struct t_mdatoms;
struct t_nrnb;
void done_lincs(Lincs* li);
/*! \brief Initialize lincs stuff */
-void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_commrec* cr, Lincs* li);
+void set_lincs(const InteractionDefinitions& idef,
+ const t_mdatoms& md,
+ bool bDynamics,
+ const t_commrec* cr,
+ Lincs* li);
/*! \brief Applies LINCS constraints.
*
#include "gromacs/mdlib/constr.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pbcutil/pbc_aiuc_cuda.cuh"
+#include "gromacs/topology/forcefieldparameters.h"
#include "gromacs/topology/ifunc.h"
#include "gromacs/topology/topology.h"
return true;
}
-void LincsGpu::set(const t_idef& idef, const t_mdatoms& md)
+void LincsGpu::set(const InteractionDefinitions& idef, const t_mdatoms& md)
{
int numAtoms = md.nr;
// List of constrained atoms (CPU memory)
std::vector<float> massFactorsHost;
// List of constrained atoms in local topology
- gmx::ArrayRef<const int> iatoms =
- constArrayRefFromArray(idef.il[F_CONSTR].iatoms, idef.il[F_CONSTR].nr);
- const int stride = NRAL(F_CONSTR) + 1;
- const int numConstraints = idef.il[F_CONSTR].nr / stride;
+ ArrayRef<const int> iatoms = idef.il[F_CONSTR].iatoms;
+ const int stride = NRAL(F_CONSTR) + 1;
+ const int numConstraints = idef.il[F_CONSTR].size() / stride;
// Early exit if no constraints
if (numConstraints == 0)
#include "gromacs/mdlib/constr.h"
#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/pbcutil/pbc_aiuc.h"
-#include "gromacs/topology/idef.h"
#include "gromacs/utility/classhelpers.h"
+class InteractionDefinitions;
+
namespace gmx
{
* \param[in] idef Local topology data to get information on constraints from.
* \param[in] md Atoms data to get atom masses from.
*/
- void set(const t_idef& idef, const t_mdatoms& md);
+ void set(const InteractionDefinitions& idef, const t_mdatoms& md);
/*! \brief
* Returns whether the maximum number of coupled constraints is supported
sfree(settled);
}
-void settle_set_constraints(settledata* settled, const t_ilist* il_settle, const t_mdatoms& mdatoms)
+void settle_set_constraints(settledata* settled, const InteractionList& il_settle, const t_mdatoms& mdatoms)
{
#if GMX_SIMD_HAVE_REAL
const int pack_size = GMX_SIMD_REAL_WIDTH;
#endif
const int nral1 = 1 + NRAL(F_SETTLE);
- int nsettle = il_settle->nr / nral1;
+ int nsettle = il_settle.size() / nral1;
settled->nsettle = nsettle;
if (nsettle > 0)
{
- const t_iatom* iatoms = il_settle->iatoms;
+ ArrayRef<const int> iatoms = il_settle.iatoms;
/* Here we initialize the normal SETTLE parameters */
if (settled->massw.mO < 0)
#include "gromacs/topology/idef.h"
-struct gmx_cmap_t;
struct gmx_mtop_t;
+struct InteractionList;
struct t_inputrec;
struct t_mdatoms;
struct t_pbc;
void settle_free(settledata* settled);
/*! \brief Set up the indices for the settle constraints */
-void settle_set_constraints(settledata* settled, const t_ilist* il_settle, const t_mdatoms& mdatoms);
+void settle_set_constraints(settledata* settled, const InteractionList& il_settle, const t_mdatoms& mdatoms);
/*! \brief Constrain coordinates using SETTLE.
* Can be called on any number of threads.
void settle_proj(settledata* settled,
ConstraintVariable econq,
int nsettle,
- const t_iatom iatoms[],
+ const int iatoms[],
const t_pbc* pbc, /* PBC data pointer, can be NULL */
ArrayRef<const RVec> x,
ArrayRef<RVec> der,
}
}
-void SettleGpu::set(const t_idef& idef, const t_mdatoms gmx_unused& md)
+void SettleGpu::set(const InteractionDefinitions& idef, const t_mdatoms gmx_unused& md)
{
- const int nral1 = 1 + NRAL(F_SETTLE);
- t_ilist il_settle = idef.il[F_SETTLE];
- t_iatom* iatoms = il_settle.iatoms;
- numSettles_ = il_settle.nr / nral1;
+ const int nral1 = 1 + NRAL(F_SETTLE);
+ const InteractionList& il_settle = idef.il[F_SETTLE];
+ ArrayRef<const int> iatoms = il_settle.iatoms;
+ numSettles_ = il_settle.size() / nral1;
reallocateDeviceBuffer(&d_atomIds_, numSettles_, &numAtomIds_, &numAtomIdsAlloc_, nullptr);
h_atomIds_.resize(numSettles_);
#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pbcutil/pbc_aiuc.h"
-#include "gromacs/topology/idef.h"
#include "gromacs/topology/topology.h"
+class InteractionDefinitions;
+
namespace gmx
{
* \param[in] idef System topology
* \param[in] md Atoms data. Can be used to update masses if needed (not used now).
*/
- void set(const t_idef& idef, const t_mdatoms& md);
+ void set(const InteractionDefinitions& idef, const t_mdatoms& md);
private:
//! GPU stream
}
}
-void make_shake_sblock_serial(shakedata* shaked, const t_idef* idef, const t_mdatoms& md)
+void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, const t_mdatoms& md)
{
int i, j, m, ncons;
int bstart, bnr;
t_blocka sblocks;
t_sortblock* sb;
- t_iatom* iatom;
int* inv_sblock;
/* Since we are processing the local topology,
* the F_CONSTRNC ilist has been concatenated to the F_CONSTR ilist.
*/
- ncons = idef->il[F_CONSTR].nr / 3;
+ ncons = idef->il[F_CONSTR].size() / 3;
init_blocka(&sblocks);
sfree(sblocks.index); // To solve memory leak
- gen_sblocks(nullptr, 0, md.homenr, idef, &sblocks, FALSE);
+ gen_sblocks(nullptr, 0, md.homenr, *idef, &sblocks, FALSE);
/*
bstart=(idef->nodeid > 0) ? blocks->multinr[idef->nodeid-1] : 0;
* sort the constraints in order of the sblock number
* and the atom numbers, really sorting a segment of the array!
*/
- iatom = idef->il[F_CONSTR].iatoms;
+ int* iatom = idef->il[F_CONSTR].iatoms.data();
snew(sb, ncons);
for (i = 0; (i < ncons); i++, iatom += 3)
{
pr_sortblock(debug, "After sorting", ncons, sb);
}
- iatom = idef->il[F_CONSTR].iatoms;
+ iatom = idef->il[F_CONSTR].iatoms.data();
for (i = 0; (i < ncons); i++, iatom += 3)
{
for (m = 0; (m < 3); m++)
}
// TODO: Check if this code is useful. It might never be called.
-void make_shake_sblock_dd(shakedata* shaked, const t_ilist* ilcon, const gmx_domdec_t* dd)
+void make_shake_sblock_dd(shakedata* shaked, const InteractionList& ilcon, const gmx_domdec_t* dd)
{
- int ncons, c, cg;
- t_iatom* iatom;
+ int ncons, c, cg;
if (dd->ncg_home + 1 > shaked->sblock_nalloc)
{
srenew(shaked->sblock, shaked->sblock_nalloc);
}
- ncons = ilcon->nr / 3;
- iatom = ilcon->iatoms;
- shaked->nblocks = 0;
- cg = 0;
+ ncons = ilcon.size() / 3;
+ const int* iatom = ilcon.iatoms.data();
+ shaked->nblocks = 0;
+ cg = 0;
for (c = 0; c < ncons; c++)
{
if (c == 0 || iatom[1] >= cg + 1)
}
//! Applies SHAKE
-static int vec_shakef(FILE* fplog,
- shakedata* shaked,
- const real invmass[],
- int ncon,
- t_iparams ip[],
- t_iatom* iatom,
- real tol,
- const rvec x[],
- rvec prime[],
- real omega,
- bool bFEP,
- real lambda,
- real scaled_lagrange_multiplier[],
- real invdt,
- rvec* v,
- bool bCalcVir,
- tensor vir_r_m_dr,
- ConstraintVariable econq)
+static int vec_shakef(FILE* fplog,
+ shakedata* shaked,
+ const real invmass[],
+ int ncon,
+ ArrayRef<const t_iparams> ip,
+ const int* iatom,
+ real tol,
+ const rvec x[],
+ rvec prime[],
+ real omega,
+ bool bFEP,
+ real lambda,
+ real scaled_lagrange_multiplier[],
+ real invdt,
+ rvec* v,
+ bool bCalcVir,
+ tensor vir_r_m_dr,
+ ConstraintVariable econq)
{
- rvec* rij;
- real * half_of_reduced_mass, *distance_squared_tolerance, *constraint_distance_squared;
- int maxnit = 1000;
- int nit = 0, ll, i, j, d, d2, type;
- t_iatom* ia;
- real L1;
- real mm = 0., tmp;
- int error = 0;
- real constraint_distance;
+ rvec* rij;
+ real *half_of_reduced_mass, *distance_squared_tolerance, *constraint_distance_squared;
+ int maxnit = 1000;
+ int nit = 0, ll, i, j, d, d2, type;
+ real L1;
+ real mm = 0., tmp;
+ int error = 0;
+ real constraint_distance;
if (ncon > shaked->nalloc)
{
distance_squared_tolerance = shaked->distance_squared_tolerance;
constraint_distance_squared = shaked->constraint_distance_squared;
- L1 = 1.0 - lambda;
- ia = iatom;
+ L1 = 1.0 - lambda;
+ const int* ia = iatom;
for (ll = 0; (ll < ncon); ll++, ia += 3)
{
type = ia[0];
}
//! Check that constraints are satisfied.
-static void check_cons(FILE* log,
- int nc,
- const rvec x[],
- rvec prime[],
- rvec v[],
- t_iparams ip[],
- t_iatom* iatom,
- const real invmass[],
- ConstraintVariable econq)
+static void check_cons(FILE* log,
+ int nc,
+ const rvec x[],
+ rvec prime[],
+ rvec v[],
+ ArrayRef<const t_iparams> ip,
+ const int* iatom,
+ const real invmass[],
+ ConstraintVariable econq)
{
- t_iatom* ia;
- int ai, aj;
- int i;
- real d, dp;
- rvec dx, dv;
+ int ai, aj;
+ int i;
+ real d, dp;
+ rvec dx, dv;
GMX_ASSERT(v, "Input has to be non-null");
fprintf(log, " i mi j mj before after should be\n");
- ia = iatom;
+ const int* ia = iatom;
for (i = 0; (i < nc); i++, ia += 3)
{
ai = ia[1];
}
//! Applies SHAKE.
-static bool bshakef(FILE* log,
- shakedata* shaked,
- const real invmass[],
- const t_idef& idef,
- const t_inputrec& ir,
- const rvec x_s[],
- rvec prime[],
- t_nrnb* nrnb,
- real lambda,
- real* dvdlambda,
- real invdt,
- rvec* v,
- bool bCalcVir,
- tensor vir_r_m_dr,
- bool bDumpOnError,
- ConstraintVariable econq)
+static bool bshakef(FILE* log,
+ shakedata* shaked,
+ const real invmass[],
+ const InteractionDefinitions& idef,
+ const t_inputrec& ir,
+ const rvec x_s[],
+ rvec prime[],
+ t_nrnb* nrnb,
+ real lambda,
+ real* dvdlambda,
+ real invdt,
+ rvec* v,
+ bool bCalcVir,
+ tensor vir_r_m_dr,
+ bool bDumpOnError,
+ ConstraintVariable econq)
{
- t_iatom* iatoms;
- real * lam, dt_2, dvdl;
- int i, n0, ncon, blen, type, ll;
- int tnit = 0, trij = 0;
+ real *lam, dt_2, dvdl;
+ int i, n0, ncon, blen, type, ll;
+ int tnit = 0, trij = 0;
- ncon = idef.il[F_CONSTR].nr / 3;
+ ncon = idef.il[F_CONSTR].size() / 3;
for (ll = 0; ll < ncon; ll++)
{
// TODO Rewrite this block so that it is obvious that i, iatoms
// and lam are all iteration variables. Is this easier if the
// sblock data structure is organized differently?
- iatoms = &(idef.il[F_CONSTR].iatoms[shaked->sblock[0]]);
- lam = shaked->scaled_lagrange_multiplier;
+ const int* iatoms = &(idef.il[F_CONSTR].iatoms[shaked->sblock[0]]);
+ lam = shaked->scaled_lagrange_multiplier;
for (i = 0; (i < shaked->nblocks);)
{
blen = (shaked->sblock[i + 1] - shaked->sblock[i]);
{
if (ir.efep != efepNO)
{
- real bondA, bondB;
+ ArrayRef<const t_iparams> iparams = idef.iparams;
+
/* TODO This should probably use invdt, so that sd integrator scaling works properly */
dt_2 = 1 / gmx::square(ir.delta_t);
dvdl = 0;
/* The vector scaled_lagrange_multiplier[ll] contains the value -2 r_ll eta_ll
(eta_ll is the estimate of the Langrangian, definition on page 336 of Ryckaert et
al 1977), so the pre-factors are already present. */
- bondA = idef.iparams[type].constr.dA;
- bondB = idef.iparams[type].constr.dB;
+ const real bondA = iparams[type].constr.dA;
+ const real bondB = iparams[type].constr.dB;
dvdl += shaked->scaled_lagrange_multiplier[ll] * dt_2 * (bondB - bondA);
}
*dvdlambda += dvdl;
return TRUE;
}
-bool constrain_shake(FILE* log,
- shakedata* shaked,
- const real invmass[],
- const t_idef& idef,
- const t_inputrec& ir,
- const rvec x_s[],
- rvec xprime[],
- rvec vprime[],
- t_nrnb* nrnb,
- real lambda,
- real* dvdlambda,
- real invdt,
- rvec* v,
- bool bCalcVir,
- tensor vir_r_m_dr,
- bool bDumpOnError,
- ConstraintVariable econq)
+bool constrain_shake(FILE* log,
+ shakedata* shaked,
+ const real invmass[],
+ const InteractionDefinitions& idef,
+ const t_inputrec& ir,
+ const rvec x_s[],
+ rvec xprime[],
+ rvec vprime[],
+ t_nrnb* nrnb,
+ real lambda,
+ real* dvdlambda,
+ real invdt,
+ rvec* v,
+ bool bCalcVir,
+ tensor vir_r_m_dr,
+ bool bDumpOnError,
+ ConstraintVariable econq)
{
if (shaked->nblocks == 0)
{
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#ifndef GMX_MDLIB_SHAKE_H
#define GMX_MDLIB_SHAKE_H
+#include "gromacs/math/vec.h"
#include "gromacs/topology/block.h"
-#include "gromacs/topology/idef.h"
+#include "gromacs/utility/real.h"
struct gmx_domdec_t;
+struct InteractionList;
+class InteractionDefinitions;
struct t_inputrec;
struct t_mdatoms;
struct t_nrnb;
void done_shake(shakedata* d);
//! Make SHAKE blocks when not using DD.
-void make_shake_sblock_serial(shakedata* shaked, const t_idef* idef, const t_mdatoms& md);
+void make_shake_sblock_serial(shakedata* shaked, InteractionDefinitions* idef, const t_mdatoms& md);
//! Make SHAKE blocks when using DD.
-void make_shake_sblock_dd(shakedata* shaked, const t_ilist* ilcon, const gmx_domdec_t* dd);
+void make_shake_sblock_dd(shakedata* shaked, const InteractionList& ilcon, const gmx_domdec_t* dd);
/*! \brief Shake all the atoms blockwise. It is assumed that all the constraints
* in the idef->shakes field are sorted, to ascending block nr. The
* sblock[n] to sblock[n+1]. Array sblock should be large enough.
* Return TRUE when OK, FALSE when shake-error
*/
-bool constrain_shake(FILE* log, /* Log file */
- shakedata* shaked, /* Total number of atoms */
- const real invmass[], /* Atomic masses */
- const t_idef& idef, /* The interaction def */
- const t_inputrec& ir, /* Input record */
- const rvec x_s[], /* Coords before update */
- rvec xprime[], /* Output coords when constraining x */
- rvec vprime[], /* Output coords when constraining v */
- t_nrnb* nrnb, /* Performance measure */
- real lambda, /* FEP lambda */
- real* dvdlambda, /* FEP force */
- real invdt, /* 1/delta_t */
- rvec* v, /* Also constrain v if v!=NULL */
- bool bCalcVir, /* Calculate r x m delta_r */
- tensor vir_r_m_dr, /* sum r x m delta_r */
- bool bDumpOnError, /* Dump debugging stuff on error*/
- ConstraintVariable econq); /* which type of constraint is occurring */
+bool constrain_shake(FILE* log, /* Log file */
+ shakedata* shaked, /* Total number of atoms */
+ const real invmass[], /* Atomic masses */
+ const InteractionDefinitions& idef, /* The interaction def */
+ const t_inputrec& ir, /* Input record */
+ const rvec x_s[], /* Coords before update */
+ rvec xprime[], /* Output coords when constraining x */
+ rvec vprime[], /* Output coords when constraining v */
+ t_nrnb* nrnb, /* Performance measure */
+ real lambda, /* FEP lambda */
+ real* dvdlambda, /* FEP force */
+ real invdt, /* 1/delta_t */
+ rvec* v, /* Also constrain v if v!=NULL */
+ bool bCalcVir, /* Calculate r x m delta_r */
+ tensor vir_r_m_dr, /* sum r x m delta_r */
+ bool bDumpOnError, /* Dump debugging stuff on error*/
+ ConstraintVariable econq); /* which type of constraint is occurring */
/*! \brief Regular iterative shake */
void cshake(const int iatom[],
*/
matrix virial = { { 0 } };
spread_vsite_f(vsite, x, fDirectVir, nullptr, stepWork.computeVirial, virial, nrnb,
- &top->idef, fr->pbcType, fr->bMolPBC, graph, box, cr, wcycle);
+ top->idef, fr->pbcType, fr->bMolPBC, graph, box, cr, wcycle);
forceWithVirial.addVirialContribution(virial);
}
/*! \brief Set up flags that have the lifetime of the domain indicating what type of work is there to compute.
*/
-static DomainLifetimeWorkload setupDomainLifetimeWorkload(const t_inputrec& inputrec,
- const t_forcerec& fr,
- const pull_t* pull_work,
- const gmx_edsam* ed,
- const t_idef& idef,
- const t_fcdata& fcd,
- const t_mdatoms& mdatoms,
+static DomainLifetimeWorkload setupDomainLifetimeWorkload(const t_inputrec& inputrec,
+ const t_forcerec& fr,
+ const pull_t* pull_work,
+ const gmx_edsam* ed,
+ const InteractionDefinitions& idef,
+ const t_fcdata& fcd,
+ const t_mdatoms& mdatoms,
const SimulationWorkload& simulationWork,
const StepWorkload& stepWork)
{
stateGpu->waitCoordinatesReadyOnHost(AtomLocality::NonLocal);
}
/* Compute the bonded and non-bonded energies and optionally forces */
- do_force_lowlevel(fr, inputrec, &(top->idef), cr, ms, nrnb, wcycle, mdatoms, x, hist, &forceOut, enerd,
+ do_force_lowlevel(fr, inputrec, top->idef, cr, ms, nrnb, wcycle, mdatoms, x, hist, &forceOut, enerd,
fcd, box, lambda.data(), graph, fr->mu_tot, stepWork, ddBalanceRegionHandler);
wallcycle_stop(wcycle, ewcFORCE);
if (vsite && !(fr->haveDirectVirialContributions && !stepWork.computeVirial))
{
rvec* fshift = as_rvec_array(forceOut.forceWithShiftForces().shiftForces().data());
- spread_vsite_f(vsite, as_rvec_array(x.unpaddedArrayRef().data()), f, fshift, FALSE, nullptr,
- nrnb, &top->idef, fr->pbcType, fr->bMolPBC, graph, box, cr, wcycle);
+ spread_vsite_f(vsite, as_rvec_array(x.unpaddedArrayRef().data()), f, fshift, FALSE,
+ nullptr, nrnb, top->idef, fr->pbcType, fr->bMolPBC, graph, box, cr, wcycle);
}
if (stepWork.computeVirial)
return nsid;
}
-void gen_sblocks(FILE* fp, int at_start, int at_end, const t_idef* idef, t_blocka* sblock, gmx_bool bSettle)
+void gen_sblocks(FILE* fp, int at_start, int at_end, const InteractionDefinitions& idef, t_blocka* sblock, gmx_bool bSettle)
{
t_graph* g;
int i, i0;
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2010,2014,2015,2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2010,2014,2015,2018,2019,2020, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#include "gromacs/utility/basedefinitions.h"
+class InteractionDefinitions;
struct t_blocka;
-struct t_idef;
-void gen_sblocks(FILE* fp, int at_start, int at_end, const t_idef* idef, t_blocka* sblock, gmx_bool bSettle);
+void gen_sblocks(FILE* fp, int at_start, int at_end, const InteractionDefinitions& idef, t_blocka* sblock, gmx_bool bSettle);
/* Generate shake blocks from the constraint list. Set bSettle to yes for shake
* blocks including settles. You normally do not want this.
*/
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
dHdLambdaRef_ = 0;
}
- // Constraints and their parameters (local topology)
- for (int i = 0; i < F_NRE; i++)
- {
- idef_.il[i].nr = 0;
- }
- idef_.il[F_CONSTR].nr = constraints.size();
-
- snew(idef_.il[F_CONSTR].iatoms, constraints.size());
int maxType = 0;
- for (index i = 0; i < ssize(constraints); i++)
+ for (index i = 0; i < ssize(constraints); i += 3)
{
- if (i % 3 == 0)
+ if (maxType < constraints.at(i))
{
- if (maxType < constraints.at(i))
- {
- maxType = constraints.at(i);
- }
+ maxType = constraints.at(i);
}
- idef_.il[F_CONSTR].iatoms[i] = constraints.at(i);
}
- snew(idef_.iparams, maxType + 1);
+ auto& iparams = mtop_.ffparams.iparams;
+ iparams.resize(maxType + 1);
for (index i = 0; i < ssize(constraints) / 3; i++)
{
- idef_.iparams[constraints.at(3 * i)].constr.dA = constraintsR0.at(constraints.at(3 * i));
- idef_.iparams[constraints.at(3 * i)].constr.dB = constraintsR0.at(constraints.at(3 * i));
+ iparams[constraints.at(3 * i)].constr.dA = constraintsR0.at(constraints.at(3 * i));
+ iparams[constraints.at(3 * i)].constr.dB = constraintsR0.at(constraints.at(3 * i));
+ }
+ idef_ = std::make_unique<InteractionDefinitions>(mtop_.ffparams);
+ for (index i = 0; i < ssize(constraints); i++)
+ {
+ idef_->il[F_CONSTR].iatoms.push_back(constraints.at(i));
}
// Constraints and their parameters (global topology)
molBlock.nmol = 1;
mtop_.molblock.push_back(molBlock);
- mtop_.natoms = numAtoms;
- mtop_.ffparams.iparams.resize(maxType + 1);
- for (int i = 0; i <= maxType; i++)
- {
- mtop_.ffparams.iparams.at(i) = idef_.iparams[i];
- }
+ mtop_.natoms = numAtoms;
mtop_.bIntermolecularInteractions = false;
// Coordinates and velocities
dHdLambda_ = 0;
}
-/*! \brief
- * Cleaning up the memory.
- */
-ConstraintsTestData::~ConstraintsTestData()
-{
- sfree(idef_.il[F_CONSTR].iatoms);
- sfree(idef_.iparams);
-}
-
} // namespace test
} // namespace gmx
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
#ifndef GMX_MDLIB_TESTS_CONSTRTESTDATA_H
#define GMX_MDLIB_TESTS_CONSTRTESTDATA_H
+#include <memory>
#include <vector>
#include "gromacs/gmxlib/nrnb.h"
//! Input record (info that usually in .mdp file)
t_inputrec ir_;
//! Local topology
- t_idef idef_;
+ std::unique_ptr<InteractionDefinitions> idef_;
//! MD atoms
t_mdatoms md_;
//! Multisim data
*
*/
void reset();
-
- /*! \brief
- * Cleaning up the memory.
- */
- ~ConstraintsTestData();
};
} // namespace test
void applyShake(ConstraintsTestData* testData, t_pbc gmx_unused pbc)
{
shakedata* shaked = shake_init();
- make_shake_sblock_serial(shaked, &testData->idef_, testData->md_);
+ make_shake_sblock_serial(shaked, testData->idef_.get(), testData->md_);
bool success = constrain_shake(
- nullptr, shaked, testData->invmass_.data(), testData->idef_, testData->ir_,
+ nullptr, shaked, testData->invmass_.data(), *testData->idef_, testData->ir_,
as_rvec_array(testData->x_.data()), as_rvec_array(testData->xPrime_.data()),
as_rvec_array(testData->xPrime2_.data()), &testData->nrnb_, testData->md_.lambda,
&testData->dHdLambda_, testData->invdt_, as_rvec_array(testData->v_.data()),
// Initialize LINCS
lincsd = init_lincs(nullptr, testData->mtop_, testData->nflexcon_, at2con_mt, false,
testData->ir_.nLincsIter, testData->ir_.nProjOrder);
- set_lincs(testData->idef_, testData->md_, EI_DYNAMICS(testData->ir_.eI), &testData->cr_, lincsd);
+ set_lincs(*testData->idef_, testData->md_, EI_DYNAMICS(testData->ir_.eI), &testData->cr_, lincsd);
// Evaluate constraints
bool success = constrain_lincs(
int numAtoms = testData->numAtoms_;
float3 *d_x, *d_xp, *d_v;
- lincsGpu->set(testData->idef_, testData->md_);
+ lincsGpu->set(*testData->idef_, testData->md_);
PbcAiuc pbcAiuc;
setPbcAiuc(pbc.ndim_ePBC, pbc.box, &pbcAiuc);
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2014,2015,2016,2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2014,2015,2016,2018,2019,2020, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
mtop_.moltype[0].atoms.atom[i * atomsPerSettle_ + 2].m = hydrogenMass_;
}
- // Reshape some data so it can be directly used by the SETTLE constraints
- ilist_ = { mtop_.moltype[0].ilist[F_SETTLE].size(), mtop_.moltype[0].ilist[F_SETTLE].iatoms.data(), 0 };
- idef_.il[F_SETTLE] = ilist_;
+ idef_ = std::make_unique<InteractionDefinitions>(mtop_.ffparams);
+ idef_->il[F_SETTLE] = mtop_.moltype[0].ilist[F_SETTLE];
}
SettleTestData::~SettleTestData()
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
gmx_mtop_t mtop_;
//! Atoms data
t_mdatoms mdatoms_;
- //! Interactions list
- t_ilist ilist_;
//! Local topology
- t_idef idef_;
+ std::unique_ptr<InteractionDefinitions> idef_;
//! Inverse timestep
const real reciprocalTimeStep_ = 1.0 / 0.002;
{
settledata* settled = settle_init(testData->mtop_);
- settle_set_constraints(settled, &testData->ilist_, testData->mdatoms_);
+ settle_set_constraints(settled, testData->idef_->il[F_SETTLE], testData->mdatoms_);
bool errorOccured;
int numThreads = 1;
auto settleGpu = std::make_unique<SettleGpu>(testData->mtop_, nullptr);
- settleGpu->set(testData->idef_, testData->mdatoms_);
+ settleGpu->set(*testData->idef_, testData->mdatoms_);
PbcAiuc pbcAiuc;
setPbcAiuc(pbc.ndim_ePBC, pbc.box, &pbcAiuc);
struct gmx_mtop_t;
enum class PbcType : int;
-struct t_idef;
+class InteractionDefinitions;
struct t_inputrec;
struct t_mdatoms;
struct t_pbc;
* \param[in] md Atoms data.
* \param[in] numTempScaleValues Number of temperature scaling groups. Zero for no temperature scaling.
*/
- void set(DeviceBuffer<RVec> d_x,
- DeviceBuffer<RVec> d_v,
- DeviceBuffer<RVec> d_f,
- const t_idef& idef,
- const t_mdatoms& md,
- int numTempScaleValues);
+ void set(DeviceBuffer<RVec> d_x,
+ DeviceBuffer<RVec> d_v,
+ DeviceBuffer<RVec> d_f,
+ const InteractionDefinitions& idef,
+ const t_mdatoms& md,
+ int numTempScaleValues);
/*! \brief
* Update PBC data.
void UpdateConstrainGpu::set(DeviceBuffer<RVec> /* d_x */,
DeviceBuffer<RVec> /* d_v */,
const DeviceBuffer<RVec> /* d_f */,
- const t_idef& /* idef */,
+ const InteractionDefinitions& /* idef */,
const t_mdatoms& /* md */,
const int /* numTempScaleValues */)
{
UpdateConstrainGpu::Impl::~Impl() {}
-void UpdateConstrainGpu::Impl::set(DeviceBuffer<RVec> d_x,
- DeviceBuffer<RVec> d_v,
- const DeviceBuffer<RVec> d_f,
- const t_idef& idef,
- const t_mdatoms& md,
- const int numTempScaleValues)
+void UpdateConstrainGpu::Impl::set(DeviceBuffer<RVec> d_x,
+ DeviceBuffer<RVec> d_v,
+ const DeviceBuffer<RVec> d_f,
+ const InteractionDefinitions& idef,
+ const t_mdatoms& md,
+ const int numTempScaleValues)
{
GMX_ASSERT(d_x != nullptr, "Coordinates device buffer should not be null.");
GMX_ASSERT(d_v != nullptr, "Velocities device buffer should not be null.");
impl_->scaleCoordinates(scalingMatrix);
}
-void UpdateConstrainGpu::set(DeviceBuffer<RVec> d_x,
- DeviceBuffer<RVec> d_v,
- const DeviceBuffer<RVec> d_f,
- const t_idef& idef,
- const t_mdatoms& md,
- const int numTempScaleValues)
+void UpdateConstrainGpu::set(DeviceBuffer<RVec> d_x,
+ DeviceBuffer<RVec> d_v,
+ const DeviceBuffer<RVec> d_f,
+ const InteractionDefinitions& idef,
+ const t_mdatoms& md,
+ const int numTempScaleValues)
{
impl_->set(d_x, d_v, d_f, idef, md, numTempScaleValues);
}
* \param[in] md Atoms data.
* \param[in] numTempScaleValues Number of temperature scaling groups. Set zero for no temperature coupling.
*/
- void set(DeviceBuffer<RVec> d_x,
- DeviceBuffer<RVec> d_v,
- const DeviceBuffer<RVec> d_f,
- const t_idef& idef,
- const t_mdatoms& md,
- const int numTempScaleValues);
+ void set(DeviceBuffer<RVec> d_x,
+ DeviceBuffer<RVec> d_v,
+ const DeviceBuffer<RVec> d_f,
+ const InteractionDefinitions& idef,
+ const t_mdatoms& md,
+ const int numTempScaleValues);
/*! \brief
* Update PBC data.
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
{
for (size_t i = 0; i < ilist.iatoms.size(); i += ilistStride(ilist))
{
- if (isConstraintFlexible(iparams.data(), ilist.iatoms[i]))
+ if (isConstraintFlexible(iparams, ilist.iatoms[i]))
{
return true;
}
}
/* Combine all constraint ilists into a single one */
- InteractionList constraintsCombined = jointConstraintList(moltype);
- t_ilist ilistsCombined[F_NRE];
- ilistsCombined[F_CONSTR].nr = constraintsCombined.iatoms.size();
- ilistsCombined[F_CONSTR].iatoms = constraintsCombined.iatoms.data();
- ilistsCombined[F_CONSTRNC].nr = 0;
+ std::array<InteractionList, F_NRE> ilistsCombined;
+ ilistsCombined[F_CONSTR] = jointConstraintList(moltype);
/* We "include" flexible constraints, but none are present (checked above) */
- const ListOfLists<int> at2con = make_at2con(moltype.atoms.nr, ilistsCombined, iparams.data(),
+ const ListOfLists<int> at2con = make_at2con(moltype.atoms.nr, ilistsCombined, iparams,
FlexibleConstraintTreatment::Include);
bool satisfiesCriteria = true;
int firstAtom = 0;
while (satisfiesCriteria && firstAtom < moltype.atoms.nr)
{
- int numAtomsInGroup = detectGroup(firstAtom, moltype, at2con, constraintsCombined);
+ int numAtomsInGroup = detectGroup(firstAtom, moltype, at2con, ilistsCombined[F_CONSTR]);
if (numAtomsInGroup == 0)
{
maxAtom = a;
}
}
- GMX_ASSERT(maxAtom >= 0 || settles.size() > 0,
+ GMX_ASSERT(maxAtom >= 0 || !settles.empty(),
"We should have at least two atoms in the group with constraints");
if (maxAtom < 0)
{
* Any remaining vsites are assigned to a separate master thread task.
*/
+using gmx::ArrayRef;
using gmx::RVec;
-static void init_ilist(t_ilist* ilist)
-{
- for (int i = 0; i < F_NRE; i++)
- {
- ilist[i].nr = 0;
- ilist[i].nalloc = 0;
- ilist[i].iatoms = nullptr;
- }
-}
-
/*! \brief List of atom indices belonging to a task */
struct AtomIndex
{
struct InterdependentTask
{
//! The interaction lists, only vsite entries are used
- t_ilist ilist[F_NRE];
+ InteractionLists ilist;
//! Thread/task-local force buffer
std::vector<RVec> force;
//! The atom indices of the vsites of our task
//! Flags if elements in force are spread to or not
std::vector<bool> use;
//! The number of entries set to true in use
- int nuse;
+ int nuse = 0;
//! Array of atoms indices, size nthreads, covering all nuse set elements in use
std::vector<AtomIndex> atomIndex;
//! List of tasks (force blocks) this task spread forces to
std::vector<int> spreadTask;
//! List of tasks that write to this tasks force block range
std::vector<int> reduceTask;
-
- InterdependentTask()
- {
- init_ilist(ilist);
- nuse = 0;
- }
};
/*! \brief Vsite thread task data structure */
//! End of atom range of this task
int rangeEnd;
//! The interaction lists, only vsite entries are used
- t_ilist ilist[F_NRE];
+ std::array<InteractionList, F_NRE> ilist;
//! Local fshift accumulation buffer
rvec fshift[SHIFTS];
//! Local virial dx*df accumulation buffer
{
rangeStart = -1;
rangeEnd = -1;
- init_ilist(ilist);
clear_rvecs(SHIFTS, fshift);
clear_mat(dxdf);
useInterdependentTask = false;
*
* \param[in] ilist The interaction list
*/
-template<typename T>
-static int vsiteIlistNrCount(const T* ilist)
+static int vsiteIlistNrCount(ArrayRef<const InteractionList> ilist)
{
int nr = 0;
for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
}
-static int constr_vsiten(const t_iatom* ia, const t_iparams ip[], rvec* x, const t_pbc* pbc)
+static int constr_vsiten(const t_iatom* ia, ArrayRef<const t_iparams> ip, rvec* x, const t_pbc* pbc)
{
rvec x1, dx;
dvec dsum;
}
}
-static void construct_vsites_thread(rvec x[],
- real dt,
- rvec* v,
- const t_iparams ip[],
- const t_ilist ilist[],
- const t_pbc* pbc_null)
+static void construct_vsites_thread(rvec x[],
+ real dt,
+ rvec* v,
+ ArrayRef<const t_iparams> ip,
+ ArrayRef<const InteractionList> ilist,
+ const t_pbc* pbc_null)
{
real inv_dt;
if (v != nullptr)
for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
{
- if (ilist[ftype].nr == 0)
+ if (ilist[ftype].empty())
{
continue;
}
{ // TODO remove me
int nra = interaction_function[ftype].nratoms;
int inc = 1 + nra;
- int nr = ilist[ftype].nr;
+ int nr = ilist[ftype].size();
- const t_iatom* ia = ilist[ftype].iatoms;
+ const t_iatom* ia = ilist[ftype].iatoms.data();
for (int i = 0; i < nr;)
{
}
}
-void construct_vsites(const gmx_vsite_t* vsite,
- rvec x[],
- real dt,
- rvec* v,
- const t_iparams ip[],
- const t_ilist ilist[],
- PbcType pbcType,
- gmx_bool bMolPBC,
- const t_commrec* cr,
- const matrix box)
+void construct_vsites(const gmx_vsite_t* vsite,
+ rvec x[],
+ real dt,
+ rvec* v,
+ ArrayRef<const t_iparams> ip,
+ ArrayRef<const InteractionList> ilist,
+ PbcType pbcType,
+ gmx_bool bMolPBC,
+ const t_commrec* cr,
+ const matrix box)
{
const bool useDomdec = (vsite != nullptr && vsite->useDomdec);
GMX_ASSERT(!useDomdec || (cr != nullptr && DOMAINDECOMP(cr)),
{
const gmx_molblock_t& molb = mtop.molblock[mb];
const gmx_moltype_t& molt = mtop.moltype[molb.type];
- if (vsiteIlistNrCount(molt.ilist.data()) > 0)
+ if (vsiteIlistNrCount(molt.ilist) > 0)
{
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.data(), ilist, PbcType::No, TRUE, nullptr, nullptr);
+ mtop.ffparams.iparams, molt.ilist, PbcType::No, TRUE, nullptr, nullptr);
atomOffset += molt.atoms.nr;
}
}
}
-static int spread_vsiten(const t_iatom ia[],
- const t_iparams ip[],
- const rvec x[],
- rvec f[],
- rvec fshift[],
- const t_pbc* pbc,
- const t_graph* g)
+static int spread_vsiten(const t_iatom ia[],
+ ArrayRef<const t_iparams> ip,
+ const rvec x[],
+ rvec f[],
+ rvec fshift[],
+ const t_pbc* pbc,
+ const t_graph* g)
{
rvec xv, dx, fi;
int n3, av, i, ai;
}
-static int vsite_count(const t_ilist* ilist, int ftype)
+static int vsite_count(ArrayRef<const InteractionList> ilist, int ftype)
{
if (ftype == F_VSITEN)
{
- return ilist[ftype].nr / 3;
+ return ilist[ftype].size() / 3;
}
else
{
- return ilist[ftype].nr / (1 + interaction_function[ftype].nratoms);
+ return ilist[ftype].size() / (1 + interaction_function[ftype].nratoms);
}
}
-static void spread_vsite_f_thread(const rvec x[],
- rvec f[],
- rvec* fshift,
- gmx_bool VirCorr,
- matrix dxdf,
- t_iparams ip[],
- const t_ilist ilist[],
- const t_graph* g,
- const t_pbc* pbc_null)
+static void spread_vsite_f_thread(const rvec x[],
+ rvec f[],
+ rvec* fshift,
+ gmx_bool VirCorr,
+ matrix dxdf,
+ ArrayRef<const t_iparams> ip,
+ ArrayRef<const InteractionList> ilist,
+ const t_graph* g,
+ const t_pbc* pbc_null)
{
const PbcMode pbcMode = getPbcMode(pbc_null);
/* We need another pbc pointer, as with charge groups we switch per vsite */
* higher type vsites from lower types */
for (int ftype = c_ftypeVsiteEnd - 1; ftype >= c_ftypeVsiteStart; ftype--)
{
- if (ilist[ftype].nr == 0)
+ if (ilist[ftype].empty())
{
continue;
}
{ // TODO remove me
int nra = interaction_function[ftype].nratoms;
int inc = 1 + nra;
- int nr = ilist[ftype].nr;
+ int nr = ilist[ftype].size();
- const t_iatom* ia = ilist[ftype].iatoms;
+ const t_iatom* ia = ilist[ftype].iatoms.data();
if (pbcMode == PbcMode::all)
{
void spread_vsite_f(const gmx_vsite_t* vsite,
const rvec* gmx_restrict x,
rvec* gmx_restrict f,
- rvec* gmx_restrict fshift,
- gmx_bool VirCorr,
- matrix vir,
- t_nrnb* nrnb,
- const t_idef* idef,
- PbcType pbcType,
- gmx_bool bMolPBC,
- const t_graph* g,
- const matrix box,
- const t_commrec* cr,
- gmx_wallcycle* wcycle)
+ rvec* gmx_restrict fshift,
+ gmx_bool VirCorr,
+ matrix vir,
+ t_nrnb* nrnb,
+ const InteractionDefinitions& idef,
+ PbcType pbcType,
+ gmx_bool bMolPBC,
+ const t_graph* g,
+ const matrix box,
+ const t_commrec* cr,
+ gmx_wallcycle* wcycle)
{
wallcycle_start(wcycle, ewcVSITESPREAD);
const bool useDomdec = vsite->useDomdec;
{
clear_mat(dxdf);
}
- spread_vsite_f_thread(x, f, fshift, VirCorr, dxdf, idef->iparams, idef->il, g, pbc_null);
+ spread_vsite_f_thread(x, f, fshift, VirCorr, dxdf, idef.iparams, idef.il, g, pbc_null);
if (VirCorr)
{
clear_mat(vsite->tData[vsite->nthreads]->dxdf);
}
spread_vsite_f_thread(x, f, fshift, VirCorr, vsite->tData[vsite->nthreads]->dxdf,
- idef->iparams, vsite->tData[vsite->nthreads]->ilist, g, pbc_null);
+ idef.iparams, vsite->tData[vsite->nthreads]->ilist, g, pbc_null);
#pragma omp parallel num_threads(vsite->nthreads)
{
copy_rvec(f[idTask->vsite[i]], idTask->force[idTask->vsite[i]]);
}
spread_vsite_f_thread(x, as_rvec_array(idTask->force.data()), fshift_t, VirCorr,
- tData.dxdf, idef->iparams, tData.idTask.ilist, g, pbc_null);
+ tData.dxdf, idef.iparams, tData.idTask.ilist, g, pbc_null);
/* We need a barrier before reducing forces below
* that have been produced by a different thread above.
}
/* Spread the vsites that spread locally only */
- spread_vsite_f_thread(x, f, fshift_t, VirCorr, tData.dxdf, idef->iparams,
+ spread_vsite_f_thread(x, f, fshift_t, VirCorr, tData.dxdf, idef.iparams,
tData.ilist, g, pbc_null);
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
dd_move_f_vsites(cr->dd, f, fshift);
}
- inc_nrnb(nrnb, eNR_VSITE2, vsite_count(idef->il, F_VSITE2));
- inc_nrnb(nrnb, eNR_VSITE2FD, vsite_count(idef->il, F_VSITE2FD));
- inc_nrnb(nrnb, eNR_VSITE3, vsite_count(idef->il, F_VSITE3));
- inc_nrnb(nrnb, eNR_VSITE3FD, vsite_count(idef->il, F_VSITE3FD));
- inc_nrnb(nrnb, eNR_VSITE3FAD, vsite_count(idef->il, F_VSITE3FAD));
- inc_nrnb(nrnb, eNR_VSITE3OUT, vsite_count(idef->il, F_VSITE3OUT));
- inc_nrnb(nrnb, eNR_VSITE4FD, vsite_count(idef->il, F_VSITE4FD));
- inc_nrnb(nrnb, eNR_VSITE4FDN, vsite_count(idef->il, F_VSITE4FDN));
- inc_nrnb(nrnb, eNR_VSITEN, vsite_count(idef->il, F_VSITEN));
+ inc_nrnb(nrnb, eNR_VSITE2, vsite_count(idef.il, F_VSITE2));
+ inc_nrnb(nrnb, eNR_VSITE2FD, vsite_count(idef.il, F_VSITE2FD));
+ inc_nrnb(nrnb, eNR_VSITE3, vsite_count(idef.il, F_VSITE3));
+ inc_nrnb(nrnb, eNR_VSITE3FD, vsite_count(idef.il, F_VSITE3FD));
+ inc_nrnb(nrnb, eNR_VSITE3FAD, vsite_count(idef.il, F_VSITE3FAD));
+ inc_nrnb(nrnb, eNR_VSITE3OUT, vsite_count(idef.il, F_VSITE3OUT));
+ inc_nrnb(nrnb, eNR_VSITE4FD, vsite_count(idef.il, F_VSITE4FD));
+ inc_nrnb(nrnb, eNR_VSITE4FDN, vsite_count(idef.il, F_VSITE4FDN));
+ inc_nrnb(nrnb, eNR_VSITEN, vsite_count(idef.il, F_VSITEN));
wallcycle_stop(wcycle, ewcVSITESPREAD);
}
* taskIndex[] is set for all vsites in our range, either to our local tasks
* or to the single last task as taskIndex[]=2*nthreads.
*/
-static void assignVsitesToThread(VsiteThread* tData,
- int thread,
- int nthread,
- int natperthread,
- gmx::ArrayRef<int> taskIndex,
- const t_ilist* ilist,
- const t_iparams* ip,
- const unsigned short* ptype)
+static void assignVsitesToThread(VsiteThread* tData,
+ int thread,
+ int nthread,
+ int natperthread,
+ gmx::ArrayRef<int> taskIndex,
+ ArrayRef<const InteractionList> ilist,
+ ArrayRef<const t_iparams> ip,
+ const unsigned short* ptype)
{
for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
{
- tData->ilist[ftype].nr = 0;
- tData->idTask.ilist[ftype].nr = 0;
+ tData->ilist[ftype].clear();
+ tData->idTask.ilist[ftype].clear();
- int nral1 = 1 + NRAL(ftype);
- int inc = nral1;
- t_iatom* iat = ilist[ftype].iatoms;
- for (int i = 0; i < ilist[ftype].nr;)
+ int nral1 = 1 + NRAL(ftype);
+ int inc = nral1;
+ const int* iat = ilist[ftype].iatoms.data();
+ for (int i = 0; i < ilist[ftype].size();)
{
if (ftype == F_VSITEN)
{
if (task == thread || task == nthread + thread)
{
/* Copy this vsite to the thread data struct of thread */
- t_ilist* il_task;
+ InteractionList* il_task;
if (task == thread)
{
il_task = &tData->ilist[ftype];
{
il_task = &tData->idTask.ilist[ftype];
}
- /* Ensure we have sufficient memory allocated */
- if (il_task->nr + inc > il_task->nalloc)
- {
- il_task->nalloc = over_alloc_large(il_task->nr + inc);
- srenew(il_task->iatoms, il_task->nalloc);
- }
/* Copy the vsite data to the thread-task local array */
- for (int j = i; j < i + inc; j++)
- {
- il_task->iatoms[il_task->nr++] = iat[j];
- }
+ il_task->push_back(iat[i], nral1 - 1, iat + i + 1);
if (task == nthread + thread)
{
/* This vsite write outside our own task force block.
}
/*! \brief Assign all vsites with taskIndex[]==task to task tData */
-static void assignVsitesToSingleTask(VsiteThread* tData,
- int task,
- gmx::ArrayRef<const int> taskIndex,
- const t_ilist* ilist,
- const t_iparams* ip)
+static void assignVsitesToSingleTask(VsiteThread* tData,
+ int task,
+ gmx::ArrayRef<const int> taskIndex,
+ ArrayRef<const InteractionList> ilist,
+ ArrayRef<const t_iparams> ip)
{
for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
{
- tData->ilist[ftype].nr = 0;
- tData->idTask.ilist[ftype].nr = 0;
+ tData->ilist[ftype].clear();
+ tData->idTask.ilist[ftype].clear();
- int nral1 = 1 + NRAL(ftype);
- int inc = nral1;
- t_iatom* iat = ilist[ftype].iatoms;
- t_ilist* il_task = &tData->ilist[ftype];
+ int nral1 = 1 + NRAL(ftype);
+ int inc = nral1;
+ const int* iat = ilist[ftype].iatoms.data();
+ InteractionList* il_task = &tData->ilist[ftype];
- for (int i = 0; i < ilist[ftype].nr;)
+ for (int i = 0; i < ilist[ftype].size();)
{
if (ftype == F_VSITEN)
{
/* Check if the vsite is assigned to our task */
if (taskIndex[iat[1 + i]] == task)
{
- /* Ensure we have sufficient memory allocated */
- if (il_task->nr + inc > il_task->nalloc)
- {
- il_task->nalloc = over_alloc_large(il_task->nr + inc);
- srenew(il_task->iatoms, il_task->nalloc);
- }
/* Copy the vsite data to the thread-task local array */
- for (int j = i; j < i + inc; j++)
- {
- il_task->iatoms[il_task->nr++] = iat[j];
- }
+ il_task->push_back(iat[i], inc - 1, iat + i + 1);
}
i += inc;
}
}
-void split_vsites_over_threads(const t_ilist* ilist, const t_iparams* ip, const t_mdatoms* mdatoms, gmx_vsite_t* vsite)
+void split_vsites_over_threads(ArrayRef<const InteractionList> ilist,
+ ArrayRef<const t_iparams> ip,
+ const t_mdatoms* mdatoms,
+ gmx_vsite_t* vsite)
{
int vsite_atom_range, natperthread;
{ // TODO remove me
if (ftype != F_VSITEN)
{
- int nral1 = 1 + NRAL(ftype);
- const t_iatom* iat = ilist[ftype].iatoms;
- for (int i = 0; i < ilist[ftype].nr; i += nral1)
+ int nral1 = 1 + NRAL(ftype);
+ ArrayRef<const int> iat = ilist[ftype].iatoms;
+ for (int i = 0; i < ilist[ftype].size(); i += nral1)
{
for (int j = i + 1; j < i + nral1; j++)
{
{
int vs_ind_end;
- const t_iatom* iat = ilist[ftype].iatoms;
+ ArrayRef<const int> iat = ilist[ftype].iatoms;
int i = 0;
- while (i < ilist[ftype].nr)
+ while (i < ilist[ftype].size())
{
/* The 3 below is from 1+NRAL(ftype)=3 */
vs_ind_end = i + ip[iat[i]].vsiten.n * 3;
for (int ftype = c_ftypeVsiteStart; ftype < c_ftypeVsiteEnd; ftype++)
{
- if (ilist[ftype].nr > 0)
+ if (!ilist[ftype].empty())
{
fprintf(debug, "%-20s thread dist:", interaction_function[ftype].longname);
for (int th = 0; th < vsite->nthreads + 1; th++)
{
- fprintf(debug, " %4d %4d ", vsite->tData[th]->ilist[ftype].nr,
- vsite->tData[th]->idTask.ilist[ftype].nr);
+ fprintf(debug, " %4d %4d ", vsite->tData[th]->ilist[ftype].size(),
+ vsite->tData[th]->idTask.ilist[ftype].size());
}
fprintf(debug, "\n");
}
struct gmx_mtop_t;
struct t_commrec;
struct t_graph;
-struct t_ilist;
+struct InteractionList;
struct t_mdatoms;
struct t_nrnb;
struct gmx_wallcycle;
* \param[in] cr The communication record
* \param[in] box The box
*/
-void construct_vsites(const gmx_vsite_t* vsite,
- rvec x[],
- real dt,
- rvec v[],
- const t_iparams ip[],
- const t_ilist ilist[],
- PbcType pbcType,
- gmx_bool bMolPBC,
- const t_commrec* cr,
- const matrix box);
+void construct_vsites(const gmx_vsite_t* vsite,
+ rvec x[],
+ real dt,
+ rvec v[],
+ gmx::ArrayRef<const t_iparams> ip,
+ gmx::ArrayRef<const InteractionList> ilist,
+ PbcType pbcType,
+ gmx_bool bMolPBC,
+ const t_commrec* cr,
+ const matrix box);
/*! \brief Create positions of vsite atoms for the whole system assuming all molecules are wholex
*
*/
void constructVsitesGlobal(const gmx_mtop_t& mtop, gmx::ArrayRef<gmx::RVec> x);
-void spread_vsite_f(const gmx_vsite_t* vsite,
- const rvec x[],
- rvec f[],
- rvec* fshift,
- gmx_bool VirCorr,
- matrix vir,
- t_nrnb* nrnb,
- const t_idef* idef,
- PbcType pbcType,
- gmx_bool bMolPBC,
- const t_graph* g,
- const matrix box,
- const t_commrec* cr,
- gmx_wallcycle* wcycle);
+void spread_vsite_f(const gmx_vsite_t* vsite,
+ const rvec x[],
+ rvec f[],
+ rvec* fshift,
+ gmx_bool VirCorr,
+ matrix vir,
+ t_nrnb* nrnb,
+ const InteractionDefinitions& idef,
+ PbcType pbcType,
+ gmx_bool bMolPBC,
+ const t_graph* g,
+ const matrix box,
+ const t_commrec* cr,
+ gmx_wallcycle* wcycle);
/* Spread the force operating on the vsite atoms on the surrounding atoms.
* If fshift!=NULL also update the shift forces.
* If VirCorr=TRUE add the virial correction for non-linear vsite constructs
*/
std::unique_ptr<gmx_vsite_t> initVsite(const gmx_mtop_t& mtop, const t_commrec* cr);
-void split_vsites_over_threads(const t_ilist* ilist,
- const t_iparams* ip,
- const t_mdatoms* mdatoms,
- gmx_vsite_t* vsite);
+void split_vsites_over_threads(gmx::ArrayRef<const InteractionList> ilist,
+ gmx::ArrayRef<const t_iparams> ip,
+ const t_mdatoms* mdatoms,
+ gmx_vsite_t* vsite);
/* Divide the vsite work-load over the threads.
* Should be called at the end of the domain decomposition.
*/
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2018,2019, by the GROMACS development team, led by
+ * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
struct SimulationGroups;
struct t_forcerec;
-struct t_idef;
struct t_inputrec;
struct t_mdatoms;
struct t_nrnb;
rvec mu_tot;
matrix pressureCouplingMu, M;
gmx_repl_ex_t repl_ex = nullptr;
- gmx_localtop_t top;
PaddedHostVector<gmx::RVec> f{};
gmx_global_stat_t gstat;
t_graph* graph = nullptr;
t_state* state;
+ gmx_localtop_t top(top_global->ffparams);
+
auto mdatoms = mdAtoms->mdatoms();
std::unique_ptr<UpdateConstrainGpu> integrator;
if (DOMAINDECOMP(cr))
{
- dd_init_local_top(*top_global, &top);
-
stateInstance = std::make_unique<t_state>();
state = stateInstance.get();
dd_init_local_state(cr->dd, state_global, state);
unsigned int force_flags;
tensor force_vir, shake_vir, total_vir, pres;
rvec mu_tot;
- gmx_localtop_t top;
PaddedHostVector<gmx::RVec> f{};
gmx_global_stat_t gstat;
t_graph* graph = nullptr;
std::unique_ptr<t_state> stateInstance;
t_state* state;
+ gmx_localtop_t top(top_global->ffparams);
+
if (DOMAINDECOMP(cr))
{
- dd_init_local_top(*top_global, &top);
-
stateInstance = std::make_unique<t_state>();
state = stateInstance.get();
dd_init_local_state(cr->dd, state_global, state);
auto mdatoms = mdAtoms->mdatoms();
if (DOMAINDECOMP(cr))
{
- top->useInDomainDecomp_ = true;
- dd_init_local_top(*top_global, top);
-
dd_init_local_state(cr->dd, state_global, &ems->s);
/* Distribute the charge groups over the nodes from the master node */
{
const char* CG = "Polak-Ribiere Conjugate Gradients";
- gmx_localtop_t top;
+ gmx_localtop_t top(top_global->ffparams);
gmx_global_stat_t gstat;
t_graph* graph;
double tmp, minstep;
{
static const char* LBFGS = "Low-Memory BFGS Minimizer";
em_state_t ems;
- gmx_localtop_t top;
+ gmx_localtop_t top(top_global->ffparams);
gmx_global_stat_t gstat;
t_graph* graph;
int ncorr, nmaxcorr, point, cp, neval, nminstep;
void LegacySimulator::do_steep()
{
const char* SD = "Steepest Descents";
- gmx_localtop_t top;
+ gmx_localtop_t top(top_global->ffparams);
gmx_global_stat_t gstat;
t_graph* graph;
real stepsize;
{
const char* NM = "Normal Mode Analysis";
int nnodes;
- gmx_localtop_t top;
+ gmx_localtop_t top(top_global->ffparams);
gmx_global_stat_t gstat;
t_graph* graph;
tensor vir, pres;
* \param[in] forceRec Force record, used for constructing vsites
* \param[in,out] graph The molecular graph, used for constructing vsites when != nullptr
*/
-static void prepareRerunState(const t_trxframe& rerunFrame,
- t_state* globalState,
- bool constructVsites,
- const gmx_vsite_t* vsite,
- const t_idef& idef,
- double timeStep,
- const t_forcerec& forceRec,
- t_graph* graph)
+static void prepareRerunState(const t_trxframe& rerunFrame,
+ t_state* globalState,
+ bool constructVsites,
+ const gmx_vsite_t* vsite,
+ const InteractionDefinitions& idef,
+ double timeStep,
+ const t_forcerec& forceRec,
+ t_graph* graph)
{
auto x = makeArrayRef(globalState->x);
auto rerunX = arrayRefFromArray(reinterpret_cast<gmx::RVec*>(rerunFrame.x), globalState->natoms);
t_trxstatus* status = nullptr;
rvec mu_tot;
t_trxframe rerun_fr;
- gmx_localtop_t top;
+ gmx_localtop_t top(top_global->ffparams);
PaddedHostVector<gmx::RVec> f{};
gmx_global_stat_t gstat;
t_graph* graph = nullptr;
if (DOMAINDECOMP(cr))
{
- dd_init_local_top(*top_global, &top);
-
stateInstance = std::make_unique<t_state>();
state = stateInstance.get();
dd_init_local_state(cr->dd, state_global, state);
ArrayRef<t_shell> shells = shfc->shells;
const int nflexcon = shfc->nflexcon;
- const t_idef* idef = &top->idef;
+ const InteractionDefinitions& idef = top->idef;
if (DOMAINDECOMP(cr))
{
if (vsite)
{
construct_vsites(vsite, as_rvec_array(pos[Min].data()), inputrec->delta_t,
- as_rvec_array(v.data()), idef->iparams, idef->il, fr->pbcType,
+ as_rvec_array(v.data()), idef.iparams, idef.il, fr->pbcType,
fr->bMolPBC, cr, box);
}
{
GMX_RELEASE_ASSERT(gmx_omp_nthreads_get(emntDefault) == 1, "TPI does not support OpenMP");
- gmx_localtop_t top;
+ gmx_localtop_t top(top_global->ffparams);
PaddedHostVector<gmx::RVec> f{};
real lambda, t, temp, beta, drmax, epot;
double embU, sum_embU, *sum_UgembU, V, V_all, VembU_all;
Constraints* constr,
gmx_vsite_t* vsite) :
globalTopology_(globalTopology),
- localTopology_(std::make_unique<gmx_localtop_t>())
+ localTopology_(std::make_unique<gmx_localtop_t>(globalTopology.ffparams))
{
- if (DOMAINDECOMP(cr))
- {
- dd_init_local_top(globalTopology, localTopology_.get());
- }
- else
+ if (!DOMAINDECOMP(cr))
{
t_graph* graph = nullptr;
// Generate and initialize new topology
mk_graph_ilist(nullptr, moltype.ilist.data(), 0, moltype.atoms.nr, FALSE, FALSE, g);
}
+t_graph* mk_graph(FILE* fplog, const InteractionDefinitions& idef, int at_start, int at_end, gmx_bool bShakeOnly, gmx_bool bSettle)
+{
+ t_graph* g;
+
+ snew(g, 1);
+
+ mk_graph_ilist(fplog, idef.il.data(), at_start, at_end, bShakeOnly, bSettle, g);
+
+ return g;
+}
+
t_graph* mk_graph(FILE* fplog, const t_idef* idef, int at_start, int at_end, gmx_bool bShakeOnly, gmx_bool bSettle)
{
t_graph* g;
struct InteractionList;
struct gmx_moltype_t;
+class InteractionDefinitions;
struct t_idef;
enum class PbcType : int;
#define SHIFT_IVEC(g, i) ((g)->ishift[i])
+t_graph* mk_graph(FILE* fplog,
+ const InteractionDefinitions& idef,
+ int at_start,
+ int at_end,
+ gmx_bool bShakeOnly,
+ gmx_bool bSettle);
+/* Build a graph from an idef description. The graph can be used
+ * to generate mol-shift indices.
+ * at_start and at_end should coincide will molecule boundaries,
+ * for the whole system this is simply 0 and natoms.
+ * If bShakeOnly, only the connections in the shake list are used.
+ * If bSettle && bShakeOnly the settles are used too.
+ */
+
t_graph* mk_graph(FILE* fplog, const struct t_idef* idef, int at_start, int at_end, gmx_bool bShakeOnly, gmx_bool bSettle);
/* Build a graph from an idef description. The graph can be used
* to generate mol-shift indices.
struct gmx_rmpbc
{
- const t_idef* idef;
- int natoms_init;
- PbcType pbcType;
- int ngraph;
- rmpbc_graph_t* graph;
+ const InteractionDefinitions* interactionDefinitions;
+ const t_idef* idef;
+ int natoms_init;
+ PbcType pbcType;
+ int ePBC;
+ int ngraph;
+ rmpbc_graph_t* graph;
};
static t_graph* gmx_rmpbc_get_graph(gmx_rmpbc_t gpbc, PbcType pbcType, int natoms)
srenew(gpbc->graph, gpbc->ngraph);
gr = &gpbc->graph[gpbc->ngraph - 1];
gr->natoms = natoms;
- gr->gr = mk_graph(nullptr, gpbc->idef, 0, natoms, FALSE, FALSE);
+ if (gpbc->interactionDefinitions)
+ {
+ gr->gr = mk_graph(nullptr, *gpbc->interactionDefinitions, 0, natoms, FALSE, FALSE);
+ }
+ else
+ {
+ gr->gr = mk_graph(nullptr, gpbc->idef, 0, natoms, FALSE, FALSE);
+ }
}
return gr->gr;
}
+gmx_rmpbc_t gmx_rmpbc_init(const InteractionDefinitions& idef, PbcType pbcType, int natoms)
+{
+ gmx_rmpbc_t gpbc;
+
+ snew(gpbc, 1);
+
+ gpbc->natoms_init = natoms;
+
+ /* This sets pbc when we now it,
+ * otherwise we guess it from the instantaneous box in the trajectory.
+ */
+ gpbc->pbcType = pbcType;
+
+ gpbc->interactionDefinitions = &idef;
+
+ return gpbc;
+}
+
gmx_rmpbc_t gmx_rmpbc_init(const t_idef* idef, PbcType pbcType, int natoms)
{
gmx_rmpbc_t gpbc;
#include "gromacs/math/vectypes.h"
+class InteractionDefinitions;
struct t_atoms;
struct t_idef;
struct t_trxframe;
typedef struct gmx_rmpbc* gmx_rmpbc_t;
+gmx_rmpbc_t gmx_rmpbc_init(const InteractionDefinitions& idef, PbcType pbcType, int natoms);
+
gmx_rmpbc_t gmx_rmpbc_init(const t_idef* idef, PbcType pbcType, int natoms);
void gmx_rmpbc_done(gmx_rmpbc_t gpbc);
}
}
-static void chk_bonds(t_idef* idef, PbcType pbcType, rvec* x, matrix box, real tol)
+static void chk_bonds(const InteractionDefinitions* idef, PbcType pbcType, rvec* x, matrix box, real tol)
{
int ftype, k, ai, aj, type;
real b0, blen, deviation;
t_pbc pbc;
rvec dx;
+ gmx::ArrayRef<const t_iparams> iparams = idef->iparams;
+
set_pbc(&pbc, pbcType, box);
for (ftype = 0; (ftype < F_NRE); ftype++)
{
if ((interaction_function[ftype].flags & IF_CHEMBOND) == IF_CHEMBOND)
{
- for (k = 0; (k < idef->il[ftype].nr);)
+ for (k = 0; (k < idef->il[ftype].size());)
{
type = idef->il[ftype].iatoms[k++];
ai = idef->il[ftype].iatoms[k++];
b0 = 0;
switch (ftype)
{
- case F_BONDS: b0 = idef->iparams[type].harmonic.rA; break;
- case F_G96BONDS: b0 = std::sqrt(idef->iparams[type].harmonic.rA); break;
- case F_MORSE: b0 = idef->iparams[type].morse.b0A; break;
- case F_CUBICBONDS: b0 = idef->iparams[type].cubic.b0; break;
- case F_CONSTR: b0 = idef->iparams[type].constr.dA; break;
+ case F_BONDS: b0 = iparams[type].harmonic.rA; break;
+ case F_G96BONDS: b0 = std::sqrt(iparams[type].harmonic.rA); break;
+ case F_MORSE: b0 = iparams[type].morse.b0A; break;
+ case F_CUBICBONDS: b0 = iparams[type].cubic.b0; break;
+ case F_CONSTR: b0 = iparams[type].constr.dA; break;
default: break;
}
if (b0 != 0)
static void chk_trj(const gmx_output_env_t* oenv, const char* fn, const char* tpr, real tol)
{
- t_trxframe fr;
- t_count count;
- t_fr_time first, last;
- int j = -1, new_natoms, natoms;
- real old_t1, old_t2;
- gmx_bool bShowTimestep = TRUE, newline = FALSE;
- t_trxstatus* status;
- gmx_mtop_t mtop;
- gmx_localtop_t top;
- t_state state;
- t_inputrec ir;
-
+ t_trxframe fr;
+ t_count count;
+ t_fr_time first, last;
+ int j = -1, new_natoms, natoms;
+ real old_t1, old_t2;
+ gmx_bool bShowTimestep = TRUE, newline = FALSE;
+ t_trxstatus* status;
+ gmx_mtop_t mtop;
+ t_state state;
+ t_inputrec ir;
+
+ std::unique_ptr<gmx_localtop_t> top;
if (tpr)
{
read_tpx_state(tpr, &ir, &state, &mtop);
- gmx_mtop_generate_local_top(mtop, &top, ir.efep != efepNO);
+ top = std::make_unique<gmx_localtop_t>(mtop.ffparams);
+ gmx_mtop_generate_local_top(mtop, top.get(), ir.efep != efepNO);
}
new_natoms = -1;
natoms = -1;
natoms = new_natoms;
if (tpr)
{
- chk_bonds(&top.idef, ir.pbcType, fr.x, fr.box, tol);
+ chk_bonds(&top->idef, ir.pbcType, fr.x, fr.box, tol);
}
if (fr.bX)
{
static void reduce_ilist(gmx::ArrayRef<const int> invindex,
const std::vector<bool>& bKeep,
- t_ilist* il,
+ InteractionList* il,
int nratoms,
const char* name)
{
- t_iatom* ia;
- int i, j, newnr;
- gmx_bool bB;
-
- if (il->nr)
+ if (!il->empty())
{
- snew(ia, il->nr);
- newnr = 0;
- for (i = 0; (i < il->nr); i += nratoms + 1)
+ std::vector<int> newAtoms(nratoms);
+ InteractionList ilReduced;
+ for (int i = 0; i < il->size(); i += nratoms + 1)
{
- bB = TRUE;
- for (j = 1; (j <= nratoms); j++)
+ bool bB = true;
+ for (int j = 0; j < nratoms; j++)
{
- bB = bB && bKeep[il->iatoms[i + j]];
+ bB = bB && bKeep[il->iatoms[i + 1 + j]];
}
if (bB)
{
- ia[newnr++] = il->iatoms[i];
- for (j = 1; (j <= nratoms); j++)
+ for (int j = 0; j < nratoms; j++)
{
- ia[newnr++] = invindex[il->iatoms[i + j]];
+ newAtoms[j] = invindex[il->iatoms[i + 1 + j]];
}
+ ilReduced.push_back(il->iatoms[i], nratoms, newAtoms.data());
}
}
- fprintf(stderr, "Reduced ilist %8s from %6d to %6d entries\n", name, il->nr / (nratoms + 1),
- newnr / (nratoms + 1));
-
- il->nr = newnr;
- for (i = 0; (i < newnr); i++)
- {
- il->iatoms[i] = ia[i];
- }
+ fprintf(stderr, "Reduced ilist %8s from %6d to %6d entries\n", name,
+ il->size() / (nratoms + 1), ilReduced.size() / (nratoms + 1));
- sfree(ia);
+ *il = std::move(ilReduced);
}
}
static void reduce_topology_x(int gnx, int index[], gmx_mtop_t* mtop, rvec x[], rvec v[])
{
- gmx_localtop_t top;
+ gmx_localtop_t top(mtop->ffparams);
gmx_mtop_generate_local_top(*mtop, &top, false);
t_atoms atoms = gmx_mtop_global_atoms(mtop);
mtop->moltype[0].excls = reduce_listoflists(invindex, bKeep, top.excls, "excls");
for (int i = 0; i < F_NRE; 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].ilist[i] = std::move(top.idef.il[i]);
}
mtop->molblock.resize(1);
#include <cstdio>
+#include "gromacs/topology/forcefieldparameters.h"
#include "gromacs/topology/ifunc.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/smalloc.h"
indent = pr_title(fp, indent, title);
pr_indent(fp, indent);
fprintf(fp, "nr: %d\n", ilist.size());
- if (ilist.size() > 0)
+ if (!ilist.empty())
{
pr_indent(fp, indent);
fprintf(fp, "iatoms:\n");
idef->iparams_fbposres = nullptr;
for (int f = 0; f < F_NRE; ++f)
{
- idef->il[f].iatoms = nullptr;
- idef->il[f].nalloc = 0;
- idef->il[f].nr = 0;
- idef->numNonperturbedInteractions[f] = 0;
+ idef->il[f].iatoms = nullptr;
+ idef->il[f].nalloc = 0;
+ idef->il[f].nr = 0;
}
- idef->cmap_grid = nullptr;
- idef->iparams_posres_nalloc = 0;
- idef->iparams_fbposres_nalloc = 0;
- idef->ilsort = ilsortUNKNOWN;
+}
+
+InteractionDefinitions::InteractionDefinitions(const gmx_ffparams_t& ffparams) :
+ iparams(ffparams.iparams),
+ functype(ffparams.functype),
+ cmap_grid(ffparams.cmap_grid)
+{
+}
+
+void InteractionDefinitions::clear()
+{
+ /* Clear the counts */
+ for (auto& ilist : il)
+ {
+ ilist.clear();
+ }
+ iparams_posres.clear();
+ iparams_fbposres.clear();
}
void done_idef(t_idef* idef)
sfree(idef->il[f].iatoms);
}
- delete idef->cmap_grid;
init_idef(idef);
}
-
-void copy_ilist(const t_ilist* src, t_ilist* dst)
-{
- dst->nr = src->nr;
- dst->nalloc = src->nr;
-
- snew(dst->iatoms, dst->nalloc);
- for (int i = 0; i < dst->nr; ++i)
- {
- dst->iatoms[i] = src->iatoms[i];
- }
-}
#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/real.h"
+struct gmx_ffparams_t;
+
typedef union t_iparams {
/* Some parameters have A and B values for free energy calculations.
* The B values are not used for regular simulations of course.
/* Returns the total number of elements in iatoms */
int size() const { return static_cast<int>(iatoms.size()); }
+ /* Returns whether the list is empty */
+ bool empty() const { return iatoms.empty(); }
+
+ /* Adds one interaction to the list */
+ template<std::size_t numAtoms>
+ void push_back(const int parameterType, const std::array<int, numAtoms>& atoms)
+ {
+ const std::size_t oldSize = iatoms.size();
+ iatoms.resize(iatoms.size() + 1 + numAtoms);
+ iatoms[oldSize] = parameterType;
+ for (std::size_t i = 0; i < numAtoms; i++)
+ {
+ iatoms[oldSize + 1 + i] = atoms[i];
+ }
+ }
+
+ /* Adds one interaction to the list */
+ void push_back(const int parameterType, const int numAtoms, const int* atoms)
+ {
+ const std::size_t oldSize = iatoms.size();
+ iatoms.resize(iatoms.size() + 1 + numAtoms);
+ iatoms[oldSize] = parameterType;
+ for (int i = 0; i < numAtoms; i++)
+ {
+ iatoms[oldSize + 1 + i] = atoms[i];
+ }
+ }
+
+ /* Appends \p ilist at the back of the list */
+ void append(const InteractionList& ilist)
+ {
+ iatoms.insert(iatoms.end(), ilist.iatoms.begin(), ilist.iatoms.end());
+ }
+
+ /* Clears the list */
+ void clear() { iatoms.clear(); }
+
/* List of interactions, see explanation further down */
std::vector<int> iatoms;
};
*
* TODO: Consider only including entries in use instead of all F_NRE
*/
-typedef std::array<InteractionList, F_NRE> InteractionLists;
+using InteractionLists = std::array<InteractionList, F_NRE>;
-/* 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.
- */
+/* Deprecated list of listed interactions */
struct t_ilist
{
/* Returns the total number of elements in iatoms */
int size() const { return nr; }
+ /* Returns whether the list is empty */
+ bool empty() const { return nr == 0; }
+
int nr;
t_iatom* iatoms;
int nalloc;
};
-/* 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.
- */
+/* TODO: Remove t_ilist and remove templating on list type in mshift.cpp */
/*
* The structs InteractionList and t_ilist defines a list of atoms with their interactions.
std::vector<InteractionListHandle> handles;
for (size_t ftype = 0; ftype < ilists.size(); ftype++)
{
- if ((interaction_function[ftype].flags & flags) && ilists[ftype].size() > 0)
+ if ((interaction_function[ftype].flags & flags) && !ilists[ftype].empty())
{
handles.push_back({ static_cast<int>(ftype), ilists[ftype].iatoms });
}
ilsortFE_SORTED
};
-typedef struct t_idef
+/* Struct with list of interaction parameters and lists of interactions
+ *
+ * TODO: Convert to a proper class with private data members so we can
+ * ensure that the free-energy sorting and sorting setting is consistent.
+ */
+class InteractionDefinitions
+{
+public:
+ /* Constructor
+ *
+ * \param[in] ffparams The interaction parameters, the lifetime of the created object should not exceed the lifetime of the passed parameters
+ */
+ InteractionDefinitions(const gmx_ffparams_t& ffparams);
+
+ // Clears data not read in from ffparams
+ void clear();
+
+ // The interaction parameters
+ const std::vector<t_iparams>& iparams;
+ // The function type per type
+ const std::vector<int>& functype;
+ // Position restraint interaction parameters
+ std::vector<t_iparams> iparams_posres;
+ // Flat-bottomed position restraint parameters
+ std::vector<t_iparams> iparams_fbposres;
+ // The list of interactions for each type. Note that some, such as LJ and COUL will have 0 entries.
+ std::array<InteractionList, F_NRE> il;
+ /* The number of non-perturbed interactions at the start of each entry in il */
+ std::array<int, F_NRE> numNonperturbedInteractions;
+ // The sorting state of interaction in il
+ int ilsort = ilsortUNKNOWN;
+ // The dihedral correction maps
+ gmx_cmap_t cmap_grid;
+};
+
+/* Deprecated interation definitions, used in t_topology */
+struct t_idef
{
int ntypes;
int atnr;
t_functype* functype;
t_iparams* iparams;
real fudgeQQ;
- gmx_cmap_t* cmap_grid;
t_iparams * iparams_posres, *iparams_fbposres;
- int iparams_posres_nalloc, iparams_fbposres_nalloc;
t_ilist il[F_NRE];
- /* The number of non-perturbed interactions at the start of each entry in il */
- int numNonperturbedInteractions[F_NRE];
- int ilsort;
-} t_idef;
+ int ilsort;
+};
/*
* The struct t_idef defines all the interactions for the complete
*/
void done_idef(t_idef* idef);
-void copy_ilist(const t_ilist* src, t_ilist* dst);
-
#endif
* The cat routines below are old code from src/kernel/topcat.c
*/
+static void ilistcat(int ftype, InteractionList* dest, const InteractionList& src, int copies, int dnum, int snum)
+{
+ int nral, c, i, a;
+
+ nral = NRAL(ftype);
+
+ size_t destIndex = dest->iatoms.size();
+ dest->iatoms.resize(dest->iatoms.size() + copies * src.size());
+
+ for (c = 0; c < copies; c++)
+ {
+ for (i = 0; i < src.size();)
+ {
+ dest->iatoms[destIndex++] = src.iatoms[i++];
+ for (a = 0; a < nral; a++)
+ {
+ dest->iatoms[destIndex++] = dnum + src.iatoms[i++];
+ }
+ }
+ dnum += snum;
+ }
+}
+
static void ilistcat(int ftype, t_ilist* dest, const InteractionList& src, int copies, int dnum, int snum)
{
int nral, c, i, a;
}
}
-static void set_posres_params(t_idef* idef, const gmx_molblock_t* molb, int i0, int a_offset)
+static const t_iparams& getIparams(const InteractionDefinitions& idef, const int index)
+{
+ return idef.iparams[index];
+}
+
+static const t_iparams& getIparams(const t_idef& idef, const int index)
+{
+ return idef.iparams[index];
+}
+
+static void resizeIParams(std::vector<t_iparams>* iparams, const int newSize)
+{
+ iparams->resize(newSize);
+}
+
+static void resizeIParams(t_iparams** iparams, const int newSize)
+{
+ srenew(*iparams, newSize);
+}
+
+template<typename IdefType>
+static void set_posres_params(IdefType* idef, const gmx_molblock_t* molb, int i0, int a_offset)
{
- t_ilist* il;
int i1, i, a_molb;
t_iparams* ip;
- il = &idef->il[F_POSRES];
- i1 = il->nr / 2;
- idef->iparams_posres_nalloc = i1;
- srenew(idef->iparams_posres, idef->iparams_posres_nalloc);
+ auto* il = &idef->il[F_POSRES];
+ i1 = il->size() / 2;
+ resizeIParams(&idef->iparams_posres, i1);
for (i = i0; i < i1; i++)
{
ip = &idef->iparams_posres[i];
/* Copy the force constants */
- *ip = idef->iparams[il->iatoms[i * 2]];
+ *ip = getIparams(*idef, il->iatoms[i * 2]);
a_molb = il->iatoms[i * 2 + 1] - a_offset;
if (molb->posres_xA.empty())
{
}
}
-static void set_fbposres_params(t_idef* idef, const gmx_molblock_t* molb, int i0, int a_offset)
+template<typename IdefType>
+static void set_fbposres_params(IdefType* idef, const gmx_molblock_t* molb, int i0, int a_offset)
{
- t_ilist* il;
int i1, i, a_molb;
t_iparams* ip;
- il = &idef->il[F_FBPOSRES];
- i1 = il->nr / 2;
- idef->iparams_fbposres_nalloc = i1;
- srenew(idef->iparams_fbposres, idef->iparams_fbposres_nalloc);
+ auto* il = &idef->il[F_FBPOSRES];
+ i1 = il->size() / 2;
+ resizeIParams(&idef->iparams_fbposres, i1);
for (i = i0; i < i1; i++)
{
ip = &idef->iparams_fbposres[i];
/* Copy the force constants */
- *ip = idef->iparams[il->iatoms[i * 2]];
+ *ip = getIparams(*idef, il->iatoms[i * 2]);
a_molb = il->iatoms[i * 2 + 1] - a_offset;
if (molb->posres_xA.empty())
{
}
}
-/*! \brief Copy idef structure from mtop.
+/*! \brief Copy parameters to idef structure from mtop.
*
- * Makes a deep copy of an idef data structure from a gmx_mtop_t.
+ * Makes a deep copy of the force field parameters data structure from a gmx_mtop_t.
* Used to initialize legacy topology types.
*
* \param[in] mtop Reference to input mtop.
* \param[in] idef Pointer to idef to populate.
- * \param[in] mergeConstr Decide if constraints will be merged.
- * \param[in] freeEnergyInteractionsAtEnd Decide if free energy stuff should
- * be added at the end.
*/
-static void copyIdefFromMtop(const gmx_mtop_t& mtop, t_idef* idef, bool freeEnergyInteractionsAtEnd, bool mergeConstr)
+static void copyFFParametersFromMtop(const gmx_mtop_t& mtop, t_idef* idef)
{
const gmx_ffparams_t* ffp = &mtop.ffparams;
{
idef->iparams = nullptr;
}
- idef->iparams_posres = nullptr;
- idef->iparams_posres_nalloc = 0;
- idef->iparams_fbposres = nullptr;
- idef->iparams_fbposres_nalloc = 0;
- idef->fudgeQQ = ffp->fudgeQQ;
- idef->cmap_grid = new gmx_cmap_t;
- *idef->cmap_grid = ffp->cmap_grid;
- idef->ilsort = ilsortUNKNOWN;
-
- for (int ftype = 0; ftype < F_NRE; ftype++)
- {
- idef->il[ftype].nr = 0;
- idef->il[ftype].nalloc = 0;
- idef->il[ftype].iatoms = nullptr;
- }
+ idef->iparams_posres = nullptr;
+ idef->iparams_fbposres = nullptr;
+ idef->fudgeQQ = ffp->fudgeQQ;
+ idef->ilsort = ilsortUNKNOWN;
+}
+/*! \brief Copy idef structure from mtop.
+ *
+ * Makes a deep copy of an idef data structure from a gmx_mtop_t.
+ * Used to initialize legacy topology types.
+ *
+ * \param[in] mtop Reference to input mtop.
+ * \param[in] idef Pointer to idef to populate.
+ * \param[in] mergeConstr Decide if constraints will be merged.
+ */
+template<typename IdefType>
+static void copyIListsFromMtop(const gmx_mtop_t& mtop, IdefType* idef, bool mergeConstr)
+{
int natoms = 0;
for (const gmx_molblock_t& molb : mtop.molblock)
{
int srcnr = molt.atoms.nr;
int destnr = natoms;
- int nposre_old = idef->il[F_POSRES].nr;
- int nfbposre_old = idef->il[F_FBPOSRES].nr;
+ int nposre_old = idef->il[F_POSRES].size();
+ int nfbposre_old = idef->il[F_FBPOSRES].size();
for (int ftype = 0; ftype < F_NRE; ftype++)
{
- if (mergeConstr && ftype == F_CONSTR && molt.ilist[F_CONSTRNC].size() > 0)
+ if (mergeConstr && ftype == F_CONSTR && !molt.ilist[F_CONSTRNC].empty())
{
/* Merge all constrains into one ilist.
* This simplifies the constraint code.
ilistcat(ftype, &idef->il[ftype], molt.ilist[ftype], molb.nmol, destnr, srcnr);
}
}
- if (idef->il[F_POSRES].nr > nposre_old)
+ if (idef->il[F_POSRES].size() > nposre_old)
{
/* Executing this line line stops gmxdump -sys working
* correctly. I'm not aware there's an elegant fix. */
set_posres_params(idef, &molb, nposre_old / 2, natoms);
}
- if (idef->il[F_FBPOSRES].nr > nfbposre_old)
+ if (idef->il[F_FBPOSRES].size() > nfbposre_old)
{
set_fbposres_params(idef, &molb, nfbposre_old / 2, natoms);
}
}
}
- if (freeEnergyInteractionsAtEnd && gmx_mtop_bondeds_free_energy(&mtop))
- {
- std::vector<real> qA(mtop.natoms);
- std::vector<real> qB(mtop.natoms);
- for (const AtomProxy atomP : AtomRange(mtop))
- {
- const t_atom& local = atomP.atom();
- int index = atomP.globalAtomNumber();
- qA[index] = local.q;
- qB[index] = local.qB;
- }
- gmx_sort_ilist_fe(idef, qA.data(), qB.data());
- }
- else
- {
- idef->ilsort = ilsortNO_FE;
- }
+ // We have not (yet) sorted free-energy interactions to the end of the ilists
+ idef->ilsort = ilsortNO_FE;
}
/*! \brief Copy atomtypes from mtop
gmx::mergeExclusions(excls, qmexcl2);
}
+static void sortFreeEnergyInteractionsAtEnd(const gmx_mtop_t& mtop, InteractionDefinitions* idef)
+{
+ std::vector<real> qA(mtop.natoms);
+ std::vector<real> qB(mtop.natoms);
+ for (const AtomProxy atomP : AtomRange(mtop))
+ {
+ const t_atom& local = atomP.atom();
+ int index = atomP.globalAtomNumber();
+ qA[index] = local.q;
+ qB[index] = local.qB;
+ }
+ gmx_sort_ilist_fe(idef, qA.data(), qB.data());
+}
+
static void gen_local_top(const gmx_mtop_t& mtop,
bool freeEnergyInteractionsAtEnd,
bool bMergeConstr,
gmx_localtop_t* top)
{
- copyAtomtypesFromMtop(mtop, &top->atomtypes);
- copyIdefFromMtop(mtop, &top->idef, freeEnergyInteractionsAtEnd, bMergeConstr);
+ copyIListsFromMtop(mtop, &top->idef, bMergeConstr);
+ if (freeEnergyInteractionsAtEnd)
+ {
+ sortFreeEnergyInteractionsAtEnd(mtop, &top->idef);
+ }
top->excls = globalExclusionLists(mtop);
if (!mtop.intermolecularExclusionGroup.empty())
{
return mols;
}
-static void gen_t_topology(const gmx_mtop_t& mtop,
- bool freeEnergyInteractionsAtEnd,
- bool bMergeConstr,
- t_topology* top)
+static void gen_t_topology(const gmx_mtop_t& mtop, bool bMergeConstr, t_topology* top)
{
copyAtomtypesFromMtop(mtop, &top->atomtypes);
- copyIdefFromMtop(mtop, &top->idef, freeEnergyInteractionsAtEnd, bMergeConstr);
+ for (int ftype = 0; ftype < F_NRE; ftype++)
+ {
+ top->idef.il[ftype].nr = 0;
+ top->idef.il[ftype].nalloc = 0;
+ top->idef.il[ftype].iatoms = nullptr;
+ }
+ copyFFParametersFromMtop(mtop, &top->idef);
+ copyIListsFromMtop(mtop, &top->idef, bMergeConstr);
top->name = mtop.name;
top->atoms = gmx_mtop_global_atoms(&mtop);
{
t_topology top;
- gen_t_topology(*mtop, false, false, &top);
+ gen_t_topology(*mtop, false, &top);
if (freeMTop)
{
}
}
-gmx_localtop_t::gmx_localtop_t()
-{
- init_idef(&idef);
- init_atomtypes(&atomtypes);
-}
-
-gmx_localtop_t::~gmx_localtop_t()
-{
- if (!useInDomainDecomp_)
- {
- done_idef(&idef);
- done_atomtypes(&atomtypes);
- }
-}
+gmx_localtop_t::gmx_localtop_t(const gmx_ffparams_t& ffparams) : idef(ffparams) {}
bool gmx_mtop_has_masses(const gmx_mtop_t* mtop)
{
struct gmx_localtop_t
{
//! Constructor used for normal operation, manages own resources.
- gmx_localtop_t();
-
- ~gmx_localtop_t();
+ gmx_localtop_t(const gmx_ffparams_t& ffparams);
//! The interaction function definition
- t_idef idef;
- //! Atomtype properties
- t_atomtypes atomtypes;
+ InteractionDefinitions idef;
//! The exclusions
gmx::ListOfLists<int> excls;
- //! Flag for domain decomposition so we don't free already freed memory.
- bool useInDomainDecomp_ = false;
};
/* The old topology struct, completely written out, used in analysis tools */
return bPert;
}
-void gmx_sort_ilist_fe(t_idef* idef, const real* qA, const real* qB)
+void gmx_sort_ilist_fe(InteractionDefinitions* idef, const real* qA, const real* qB)
{
int ftype, nral, i, ic, ib, a;
- t_ilist* ilist;
- t_iatom* iatoms;
t_iatom* iabuf;
int iabuf_nalloc;
iabuf_nalloc = 0;
iabuf = nullptr;
- const t_iparams* iparams = idef->iparams;
-
for (ftype = 0; ftype < F_NRE; ftype++)
{
if (interaction_function[ftype].flags & IF_BOND)
{
- ilist = &idef->il[ftype];
- iatoms = ilist->iatoms;
- nral = NRAL(ftype);
- ic = 0;
- ib = 0;
- i = 0;
- while (i < ilist->nr)
+ InteractionList* ilist = &idef->il[ftype];
+ int* iatoms = ilist->iatoms.data();
+ nral = NRAL(ftype);
+ ic = 0;
+ ib = 0;
+ i = 0;
+ while (i < ilist->size())
{
/* Check if this interaction is perturbed */
- if (ip_q_pert(ftype, iatoms + i, iparams, qA, qB))
+ if (ip_q_pert(ftype, iatoms + i, idef->iparams.data(), qA, qB))
{
/* Copy to the perturbed buffer */
if (ib + 1 + nral > iabuf_nalloc)
}
}
}
- /* Now we now the number of non-perturbed interactions */
+ /* Now we know the number of non-perturbed interactions */
idef->numNonperturbedInteractions[ftype] = ic;
/* Copy the buffer with perturbed interactions to the ilist */
{
const int numNonperturbed = idef->numNonperturbedInteractions[ftype];
fprintf(debug, "%s non-pert %d pert %d\n", interaction_function[ftype].longname,
- numNonperturbed, ilist->nr - numNonperturbed);
+ numNonperturbed, ilist->size() - numNonperturbed);
}
}
}
#include "gromacs/utility/real.h"
struct gmx_mtop_t;
-struct t_idef;
+class InteractionDefinitions;
/* Returns if there are perturbed bonded interactions */
gmx_bool gmx_mtop_bondeds_free_energy(const struct gmx_mtop_t* mtop);
/* Sort all the bonded ilists in idef to have the perturbed ones at the end
* and set nr_nr_nonperturbed in ilist.
*/
-void gmx_sort_ilist_fe(struct t_idef* idef, const real* qA, const real* qB);
+void gmx_sort_ilist_fe(InteractionDefinitions* idef, const real* qA, const real* qB);
#endif
// Do lazy initialization
if (expandedTopology_ == nullptr && hasTopology())
{
- expandedTopology_ = std::make_unique<gmx_localtop_t>();
+ expandedTopology_ = std::make_unique<gmx_localtop_t>(mtop_->ffparams);
gmx_mtop_generate_local_top(*mtop_, expandedTopology_.get(), false);
}
{
GMX_RELEASE_ASSERT(topInfo.hasTopology(), "Cannot remove PBC without a topology");
- return gmx_rmpbc_init(&topInfo.expandedTopology()->idef, topInfo.pbcType(), topInfo.mtop()->natoms);
+ return gmx_rmpbc_init(topInfo.expandedTopology()->idef, topInfo.pbcType(), topInfo.mtop()->natoms);
}
} // namespace gmx
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff