#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/gmxassert.h"
-#include "gromacs/utility/smalloc.h"
+#include "gromacs/utility/listoflists.h"
#include "domdec_internal.h"
#include "domdec_specatomcomm.h"
+using gmx::ListOfLists;
+
/*! \brief Struct used during constraint setup with domain decomposition */
struct gmx_domdec_constraints_t
{
int nrec,
gmx::ArrayRef<const int> ia1,
gmx::ArrayRef<const int> ia2,
- const t_blocka* at2con,
+ const ListOfLists<int>& at2con,
const gmx_ga2la_t& ga2la,
gmx_bool bHomeConnect,
gmx_domdec_constraints_t* dc,
t_ilist* il_local,
std::vector<int>* ireq)
{
- int a1_gl, a2_gl, i, coni, b;
- const t_iatom* iap;
-
if (!dc->gc_req[con_offset + con])
{
/* Add this non-home constraint to the list */
il_local->nalloc = over_alloc_dd(il_local->nr + 3);
srenew(il_local->iatoms, il_local->nalloc);
}
- iap = constr_iatomptr(ia1, ia2, con);
+ const int* iap = constr_iatomptr(ia1, ia2, con);
il_local->iatoms[il_local->nr++] = iap[0];
- a1_gl = offset + iap[1];
- a2_gl = offset + iap[2];
+ const int a1_gl = offset + iap[1];
+ const int a2_gl = offset + iap[2];
/* The following indexing code can probably be optizimed */
if (const int* a_loc = ga2la.findHome(a1_gl))
{
if (nrec > 0)
{
- for (i = at2con->index[a]; i < at2con->index[a + 1]; i++)
+ /* Loop over the constraint connected to atom a */
+ for (const int coni : at2con[a])
{
- coni = at2con->a[i];
if (coni != con)
{
/* Walk further */
- iap = constr_iatomptr(ia1, ia2, coni);
+ const int* iap = constr_iatomptr(ia1, ia2, coni);
+ int b;
if (a == iap[1])
{
b = iap[2];
}
/*! \brief Looks up constraint for the local atoms */
-static void atoms_to_constraints(gmx_domdec_t* dd,
- const gmx_mtop_t* mtop,
- const int* cginfo,
- gmx::ArrayRef<const t_blocka> at2con_mt,
- int nrec,
- t_ilist* ilc_local,
- std::vector<int>* ireq)
+static void atoms_to_constraints(gmx_domdec_t* dd,
+ const gmx_mtop_t* mtop,
+ const int* cginfo,
+ gmx::ArrayRef<const ListOfLists<int>> at2con_mt,
+ int nrec,
+ t_ilist* ilc_local,
+ std::vector<int>* ireq)
{
- const t_blocka* at2con;
- int b_lo, offset, b_mol, i, con, con_offset;
-
gmx_domdec_constraints_t* dc = dd->constraints;
gmx_domdec_specat_comm_t* dcc = dd->constraint_comm;
* This is only required for the global index to make sure
* that we use each constraint only once.
*/
- con_offset = dc->molb_con_offset[mb] + molnr * dc->molb_ncon_mol[mb];
+ const int con_offset = dc->molb_con_offset[mb] + molnr * dc->molb_ncon_mol[mb];
/* The global atom number offset for this molecule */
- offset = a_gl - a_mol;
- at2con = &at2con_mt[molb.type];
- for (i = at2con->index[a_mol]; i < at2con->index[a_mol + 1]; i++)
+ const int offset = a_gl - a_mol;
+ /* Loop over the constraints connected to atom a_mol in the molecule */
+ const auto& at2con = at2con_mt[molb.type];
+ for (const int con : at2con[a_mol])
{
- con = at2con->a[i];
const int* iap = constr_iatomptr(ia1, ia2, con);
+ int b_mol;
if (a_mol == iap[1])
{
b_mol = iap[2];
ilc_local->nalloc = over_alloc_dd(ilc_local->nr + 3);
srenew(ilc_local->iatoms, ilc_local->nalloc);
}
- b_lo = *a_loc;
+ 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);
int nrec,
t_ilist* il_local)
{
- gmx_domdec_constraints_t* dc;
- t_ilist * ilc_local, *ils_local;
- std::vector<int>* ireq;
- gmx::ArrayRef<const t_blocka> at2con_mt;
- gmx::HashedMap<int>* ga2la_specat;
- int at_end, i, j;
- t_iatom* iap;
+ gmx_domdec_constraints_t* dc;
+ t_ilist * 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
dc->ncon = 0;
ilc_local->nr = 0;
+ gmx::ArrayRef<const ListOfLists<int>> at2con_mt;
+ std::vector<int>* ireq = nullptr;
if (dd->constraint_comm)
{
// TODO Perhaps gmx_domdec_constraints_t should keep a valid constr?
ireq = &dc->requestedGlobalAtomIndices[0];
ireq->clear();
}
- else
- {
- // Currently unreachable
- at2con_mt = {};
- ireq = nullptr;
- }
gmx::ArrayRef<const std::vector<int>> at2settle_mt;
/* When settle works inside charge groups, we assigned them already */
#include "gromacs/utility/futil.h"
#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/keyvaluetreebuilder.h"
+#include "gromacs/utility/listoflists.h"
#include "gromacs/utility/mdmodulenotification.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/snprintf.h"
const t_atom* atom = molt.atoms.atom;
- t_blocka atomToConstraints =
+ const auto atomToConstraints =
gmx::make_at2con(molt, iparams, gmx::FlexibleConstraintTreatment::Exclude);
bool haveDecoupledMode = false;
int a1 = il.iatoms[1 + i + 1];
int a2 = il.iatoms[1 + i + 2];
if ((atom[a0].m > atom[a2].m * massFactorThreshold || atom[a2].m > atom[a0].m * massFactorThreshold)
- && atomToConstraints.index[a0 + 1] - atomToConstraints.index[a0] == 1
- && atomToConstraints.index[a2 + 1] - atomToConstraints.index[a2] == 1
- && atomToConstraints.index[a1 + 1] - atomToConstraints.index[a1] >= 3)
+ && atomToConstraints[a0].ssize() == 1 && atomToConstraints[a2].ssize() == 1
+ && atomToConstraints[a1].ssize() >= 3)
{
- int constraint0 = atomToConstraints.a[atomToConstraints.index[a0]];
- int constraint2 = atomToConstraints.a[atomToConstraints.index[a2]];
+ int constraint0 = atomToConstraints[a0][0];
+ int constraint2 = atomToConstraints[a2][0];
bool foundAtom0 = false;
bool foundAtom2 = false;
- for (int conIndex = atomToConstraints.index[a1];
- conIndex < atomToConstraints.index[a1 + 1]; conIndex++)
+ for (const int constraint : atomToConstraints[a1])
{
- if (atomToConstraints.a[conIndex] == constraint0)
+ if (constraint == constraint0)
{
foundAtom0 = true;
}
- if (atomToConstraints.a[conIndex] == constraint2)
+ if (constraint == constraint2)
{
foundAtom2 = true;
}
}
}
- done_blocka(&atomToConstraints);
-
return haveDecoupledMode;
}
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/pulling/pull.h"
#include "gromacs/timing/wallcycle.h"
-#include "gromacs/topology/block.h"
#include "gromacs/topology/ifunc.h"
#include "gromacs/topology/mtop_lookup.h"
#include "gromacs/topology/mtop_util.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/listoflists.h"
#include "gromacs/utility/pleasecite.h"
#include "gromacs/utility/smalloc.h"
#include "gromacs/utility/txtdump.h"
//! The number of flexible constraints.
int nflexcon = 0;
//! A list of atoms to constraints for each moleculetype.
- std::vector<t_blocka> at2con_mt;
+ std::vector<ListOfLists<int>> at2con_mt;
//! A list of atoms to settles for each moleculetype
std::vector<std::vector<int>> at2settle_mt;
//! LINCS data.
* \param[in] numAtoms The number of atoms to construct the list for
* \param[in] ilists The interaction lists, size F_NRE
* \param[in] iparams Interaction parameters, can be null when
- * flexibleConstraintTreatment=Include \param[in] flexibleConstraintTreatment The flexible
- * constraint treatment, see enum above \returns a block struct with all constraints for each atom
+ * \p flexibleConstraintTreatment==Include
+ * \param[in] flexibleConstraintTreatment The flexible constraint treatment,
+ * see enum above
+ *
+ * \returns a block struct with all constraints for each atom
*/
template<typename T>
-static t_blocka makeAtomsToConstraintsList(int numAtoms,
- const T* ilists,
- const t_iparams* iparams,
- FlexibleConstraintTreatment flexibleConstraintTreatment)
+static ListOfLists<int> makeAtomsToConstraintsList(int numAtoms,
+ const T* ilists,
+ const t_iparams* iparams,
+ FlexibleConstraintTreatment flexibleConstraintTreatment)
{
GMX_ASSERT(flexibleConstraintTreatment == FlexibleConstraintTreatment::Include || iparams != nullptr,
"With flexible constraint detection we need valid iparams");
}
}
- t_blocka at2con;
- at2con.nr = numAtoms;
- at2con.nalloc_index = at2con.nr + 1;
- snew(at2con.index, at2con.nalloc_index);
- at2con.index[0] = 0;
+ std::vector<int> listRanges(numAtoms + 1);
for (int a = 0; a < numAtoms; a++)
{
- at2con.index[a + 1] = at2con.index[a] + count[a];
- count[a] = 0;
+ listRanges[a + 1] = listRanges[a] + count[a];
+ count[a] = 0;
}
- at2con.nra = at2con.index[at2con.nr];
- at2con.nalloc_a = at2con.nra;
- snew(at2con.a, at2con.nalloc_a);
+ std::vector<int> elements(listRanges[numAtoms]);
/* The F_CONSTRNC constraints have constraint numbers
* that continue after the last F_CONSTR constraint.
{
for (int j = 1; j < 3; j++)
{
- int a = ilist.iatoms[i + j];
- at2con.a[at2con.index[a] + count[a]++] = numConstraints;
+ const int a = ilist.iatoms[i + j];
+ elements[listRanges[a] + count[a]++] = numConstraints;
}
}
numConstraints++;
}
}
- return at2con;
+ return ListOfLists<int>(std::move(listRanges), std::move(elements));
}
-t_blocka make_at2con(int numAtoms,
- const t_ilist* ilist,
- const t_iparams* iparams,
- FlexibleConstraintTreatment flexibleConstraintTreatment)
+ListOfLists<int> make_at2con(int numAtoms,
+ const t_ilist* ilist,
+ const t_iparams* iparams,
+ FlexibleConstraintTreatment flexibleConstraintTreatment)
{
return makeAtomsToConstraintsList(numAtoms, ilist, iparams, flexibleConstraintTreatment);
}
-t_blocka make_at2con(const gmx_moltype_t& moltype,
- gmx::ArrayRef<const t_iparams> iparams,
- FlexibleConstraintTreatment flexibleConstraintTreatment)
+ListOfLists<int> make_at2con(const gmx_moltype_t& moltype,
+ gmx::ArrayRef<const t_iparams> iparams,
+ FlexibleConstraintTreatment flexibleConstraintTreatment)
{
return makeAtomsToConstraintsList(moltype.atoms.nr, moltype.ilist.data(), iparams.data(),
flexibleConstraintTreatment);
* indices to constraint indices.
*
* Note that flexible constraints are only enabled with a dynamical integrator. */
-static std::vector<t_blocka> makeAtomToConstraintMappings(const gmx_mtop_t& mtop,
- FlexibleConstraintTreatment flexibleConstraintTreatment)
+static std::vector<ListOfLists<int>> makeAtomToConstraintMappings(const gmx_mtop_t& mtop,
+ FlexibleConstraintTreatment flexibleConstraintTreatment)
{
- std::vector<t_blocka> mapping;
+ std::vector<ListOfLists<int>> mapping;
mapping.reserve(mtop.moltype.size());
for (const gmx_moltype_t& moltype : mtop.moltype)
{
Constraints::Impl::~Impl()
{
- for (auto blocka : at2con_mt)
- {
- done_blocka(&blocka);
- }
if (bSettleErrorHasOccurred != nullptr)
{
sfree(bSettleErrorHasOccurred);
impl_->ed = ed;
}
-ArrayRef<const t_blocka> Constraints::atom2constraints_moltype() const
+ArrayRef<const ListOfLists<int>> Constraints::atom2constraints_moltype() const
{
return impl_->at2con_mt;
}
struct gmx_multisim_t;
struct gmx_wallcycle;
struct pull_t;
-struct t_blocka;
struct t_commrec;
struct t_ilist;
struct t_inputrec;
{
template<typename T>
class ArrayRefWithPadding;
+template<typename>
+class ListOfLists;
//! Describes supported flavours of constrained updates.
enum class ConstraintVariable : int
//! Links the essentialdynamics and constraint code.
void saveEdsamPointer(gmx_edsam* ed);
//! Getter for use by domain decomposition.
- ArrayRef<const t_blocka> atom2constraints_moltype() const;
+ ArrayRef<const ListOfLists<int>> atom2constraints_moltype() const;
//! Getter for use by domain decomposition.
ArrayRef<const std::vector<int>> atom2settle_moltype() const;
/*! \brief Returns the flexible constraint treatment depending on whether the integrator is dynamic */
FlexibleConstraintTreatment flexibleConstraintTreatment(bool haveDynamicsIntegrator);
-/*! \brief Returns a block struct to go from atoms to constraints
+/*! \brief Returns a ListOfLists object to go from atoms to constraints
*
- * The block struct will contain constraint indices with lower indices
+ * The object will contain constraint indices with lower indices
* directly matching the order in F_CONSTR and higher indices matching
* the order in F_CONSTRNC offset by the number of constraints in F_CONSTR.
*
* \param[in] moltype The molecule data
* \param[in] iparams Interaction parameters, can be null when
- * flexibleConstraintTreatment=Include \param[in] flexibleConstraintTreatment The flexible
- * constraint treatment, see enum above \returns a block struct with all constraints for each atom
+ * \p flexibleConstraintTreatment==Include
+ * \param[in] flexibleConstraintTreatment The flexible constraint treatment,
+ * see enum above
+ *
+ * \returns a ListOfLists object with all constraints for each atom
*/
-t_blocka make_at2con(const gmx_moltype_t& moltype,
- gmx::ArrayRef<const t_iparams> iparams,
- FlexibleConstraintTreatment flexibleConstraintTreatment);
+ListOfLists<int> make_at2con(const gmx_moltype_t& moltype,
+ gmx::ArrayRef<const t_iparams> iparams,
+ FlexibleConstraintTreatment flexibleConstraintTreatment);
-/*! \brief Returns a block struct to go from atoms to constraints
+/*! \brief Returns a ListOfLists object to go from atoms to constraints
*
- * The block struct will contain constraint indices with lower indices
+ * The object will contain constraint indices with lower indices
* directly matching the order in F_CONSTR and higher indices matching
* the order in F_CONSTRNC offset by the number of constraints in F_CONSTR.
*
* \param[in] numAtoms The number of atoms to construct the list for
* \param[in] ilist Interaction list, size F_NRE
* \param[in] iparams Interaction parameters, can be null when
- * flexibleConstraintTreatment=Include \param[in] flexibleConstraintTreatment The flexible
- * constraint treatment, see enum above \returns a block struct with all constraints for each atom
+ * \p flexibleConstraintTreatment==Include
+ * \param[in] flexibleConstraintTreatment The flexible constraint treatment,
+ * see enum above
+ *
+ * \returns a ListOfLists object with all constraints for each atom
*/
-t_blocka make_at2con(int numAtoms,
- const t_ilist* ilist,
- const t_iparams* iparams,
- FlexibleConstraintTreatment flexibleConstraintTreatment);
-
-/*! \brief Returns an array of atom to constraints lists for the moltypes */
-const t_blocka* atom2constraints_moltype(const Constraints* constr);
+ListOfLists<int> make_at2con(int numAtoms,
+ const t_ilist* ilist,
+ 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);
#include "gromacs/mdlib/constr.h"
#include "gromacs/mdtypes/inputrec.h"
-#include "gromacs/topology/block.h"
#include "gromacs/topology/mtop_util.h"
#include "gromacs/utility/basedefinitions.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/listoflists.h"
#include "gromacs/utility/logger.h"
#include "gromacs/utility/real.h"
-#include "gromacs/utility/smalloc.h"
namespace gmx
{
//! Recursing function to help find all adjacent constraints.
-static void constr_recur(const t_blocka* at2con,
+static void constr_recur(const ListOfLists<int>& at2con,
const InteractionLists& ilist,
gmx::ArrayRef<const t_iparams> iparams,
gmx_bool bTopB,
int at,
int depth,
int nc,
- int* path,
+ ArrayRef<int> path,
real r0,
real r1,
real* r2max,
int* count)
{
- int c, con, a1;
gmx_bool bUse;
real len, rn0, rn1;
gmx::ArrayRef<const int> ia2 = ilist[F_CONSTRNC].iatoms;
/* Loop over all constraints connected to this atom */
- for (c = at2con->index[at]; c < at2con->index[at + 1]; c++)
+ for (const int con : at2con[at])
{
- con = at2con->a[c];
/* Do not walk over already used constraints */
bUse = TRUE;
- for (a1 = 0; a1 < depth; a1++)
+ for (int a1 = 0; a1 < depth; a1++)
{
if (con == path[a1])
{
fprintf(debug,
"Found longer constraint distance: r0 %5.3f r1 %5.3f rmax %5.3f\n", rn0,
rn1, sqrt(*r2max));
- for (a1 = 0; a1 < depth; a1++)
+ for (int a1 = 0; a1 < depth; a1++)
{
fprintf(debug, " %d %5.3f", path[a1],
iparams[constr_iatomptr(ia1, ia2, con)[0]].constr.dA);
*/
if (depth + 1 < nc && *count < 1000 * nc)
{
+ int a1;
if (ia[1] == at)
{
a1 = ia[2];
gmx::ArrayRef<const t_iparams> iparams,
const t_inputrec* ir)
{
- int natoms, *path, at, count;
+ int natoms, at, count;
- t_blocka at2con;
- real r0, r1, r2maxA, r2maxB, rmax, lam0, lam1;
+ real r0, r1, r2maxA, r2maxB, rmax, lam0, lam1;
if (molt->ilist[F_CONSTR].size() == 0 && molt->ilist[F_CONSTRNC].size() == 0)
{
natoms = molt->atoms.nr;
- at2con = make_at2con(*molt, iparams, flexibleConstraintTreatment(EI_DYNAMICS(ir->eI)));
- snew(path, 1 + ir->nProjOrder);
+ const ListOfLists<int> at2con =
+ make_at2con(*molt, iparams, flexibleConstraintTreatment(EI_DYNAMICS(ir->eI)));
+ std::vector<int> path(1 + ir->nProjOrder);
for (at = 0; at < 1 + ir->nProjOrder; at++)
{
path[at] = -1;
r1 = 0;
count = 0;
- constr_recur(&at2con, molt->ilist, iparams, FALSE, at, 0, 1 + ir->nProjOrder, path, r0, r1,
+ constr_recur(at2con, molt->ilist, iparams, FALSE, at, 0, 1 + ir->nProjOrder, path, r0, r1,
&r2maxA, &count);
}
if (ir->efep == efepNO)
r0 = 0;
r1 = 0;
count = 0;
- constr_recur(&at2con, molt->ilist, iparams, TRUE, at, 0, 1 + ir->nProjOrder, path, r0,
+ constr_recur(at2con, molt->ilist, iparams, TRUE, at, 0, 1 + ir->nProjOrder, path, r0,
r1, &r2maxB, &count);
}
lam0 = ir->fepvals->init_lambda;
}
}
- done_blocka(&at2con);
- sfree(path);
-
return rmax;
}
#include "gromacs/simd/simd.h"
#include "gromacs/simd/simd_math.h"
#include "gromacs/simd/vector_operations.h"
-#include "gromacs/topology/block.h"
#include "gromacs/topology/mtop_util.h"
#include "gromacs/utility/alignedallocator.h"
#include "gromacs/utility/arrayref.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
#include "gromacs/utility/gmxomp.h"
+#include "gromacs/utility/listoflists.h"
#include "gromacs/utility/pleasecite.h"
using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
}
//! Finds all triangles of atoms that share constraints to a central atom.
-static int count_triangle_constraints(const InteractionLists& ilist, const t_blocka& at2con)
+static int count_triangle_constraints(const InteractionLists& ilist, const ListOfLists<int>& at2con)
{
- int ncon1, ncon_tot;
- int c0, n1, c1, ac1, n2, c2;
- int ncon_triangle;
-
- ncon1 = ilist[F_CONSTR].size() / 3;
- ncon_tot = ncon1 + ilist[F_CONSTRNC].size() / 3;
+ const int ncon1 = ilist[F_CONSTR].size() / 3;
+ const int ncon_tot = ncon1 + ilist[F_CONSTRNC].size() / 3;
gmx::ArrayRef<const int> ia1 = ilist[F_CONSTR].iatoms;
gmx::ArrayRef<const int> ia2 = ilist[F_CONSTRNC].iatoms;
- ncon_triangle = 0;
- for (c0 = 0; c0 < ncon_tot; c0++)
+ int ncon_triangle = 0;
+ for (int c0 = 0; c0 < ncon_tot; c0++)
{
bool bTriangle = FALSE;
const int* iap = constr_iatomptr(ia1, ia2, c0);
const int a00 = iap[1];
const int a01 = iap[2];
- for (n1 = at2con.index[a01]; n1 < at2con.index[a01 + 1]; n1++)
+ for (const int c1 : at2con[a01])
{
- c1 = at2con.a[n1];
if (c1 != c0)
{
const int* iap = constr_iatomptr(ia1, ia2, c1);
const int a10 = iap[1];
const int a11 = iap[2];
+ int ac1;
if (a10 == a01)
{
ac1 = a11;
{
ac1 = a10;
}
- for (n2 = at2con.index[ac1]; n2 < at2con.index[ac1 + 1]; n2++)
+ for (const int c2 : at2con[ac1])
{
- c2 = at2con.a[n2];
if (c2 != c0 && c2 != c1)
{
const int* iap = constr_iatomptr(ia1, ia2, c2);
}
//! Finds sequences of sequential constraints.
-static bool more_than_two_sequential_constraints(const InteractionLists& ilist, const t_blocka& at2con)
+static bool more_than_two_sequential_constraints(const InteractionLists& ilist, const ListOfLists<int>& at2con)
{
- int ncon1, ncon_tot, c;
- bool bMoreThanTwoSequentialConstraints;
-
- ncon1 = ilist[F_CONSTR].size() / 3;
- ncon_tot = ncon1 + ilist[F_CONSTRNC].size() / 3;
+ const int ncon1 = ilist[F_CONSTR].size() / 3;
+ const int ncon_tot = ncon1 + ilist[F_CONSTRNC].size() / 3;
gmx::ArrayRef<const int> ia1 = ilist[F_CONSTR].iatoms;
gmx::ArrayRef<const int> ia2 = ilist[F_CONSTRNC].iatoms;
- bMoreThanTwoSequentialConstraints = FALSE;
- for (c = 0; c < ncon_tot && !bMoreThanTwoSequentialConstraints; c++)
+ for (int c = 0; c < ncon_tot; c++)
{
const int* iap = constr_iatomptr(ia1, ia2, c);
const int a1 = iap[1];
const int a2 = iap[2];
/* Check if this constraint has constraints connected at both atoms */
- if (at2con.index[a1 + 1] - at2con.index[a1] > 1 && at2con.index[a2 + 1] - at2con.index[a2] > 1)
+ if (at2con[a1].ssize() > 1 && at2con[a2].ssize() > 1)
{
- bMoreThanTwoSequentialConstraints = TRUE;
+ return true;
}
}
- return bMoreThanTwoSequentialConstraints;
+ return false;
}
-Lincs* init_lincs(FILE* fplog,
- const gmx_mtop_t& mtop,
- int nflexcon_global,
- ArrayRef<const t_blocka> at2con,
- bool bPLINCS,
- int nIter,
- int nProjOrder)
+Lincs* init_lincs(FILE* fplog,
+ const gmx_mtop_t& mtop,
+ int nflexcon_global,
+ ArrayRef<const ListOfLists<int>> atomToConstraintsPerMolType,
+ bool bPLINCS,
+ int nIter,
+ int nProjOrder)
{
// TODO this should become a unique_ptr
Lincs* li;
li->max_connect = 0;
for (size_t mt = 0; mt < mtop.moltype.size(); mt++)
{
+ const auto& at2con = atomToConstraintsPerMolType[mt];
for (int a = 0; a < mtop.moltype[mt].atoms.nr; a++)
{
- li->max_connect = std::max(li->max_connect, at2con[mt].index[a + 1] - at2con[mt].index[a]);
+ li->max_connect = std::max(li->max_connect, int(at2con[a].ssize()));
}
}
bMoreThanTwoSeq = FALSE;
for (const gmx_molblock_t& molb : mtop.molblock)
{
- const gmx_moltype_t& molt = mtop.moltype[molb.type];
+ const gmx_moltype_t& molt = mtop.moltype[molb.type];
+ const auto& at2con = atomToConstraintsPerMolType[molb.type];
- li->ncg_triangle += molb.nmol * count_triangle_constraints(molt.ilist, at2con[molb.type]);
+ li->ncg_triangle += molb.nmol * count_triangle_constraints(molt.ilist, at2con);
- if (!bMoreThanTwoSeq && more_than_two_sequential_constraints(molt.ilist, at2con[molb.type]))
+ if (!bMoreThanTwoSeq && more_than_two_sequential_constraints(molt.ilist, at2con))
{
bMoreThanTwoSeq = TRUE;
}
}
//! Assign a constraint.
-static void assign_constraint(Lincs* li, int constraint_index, int a1, int a2, real lenA, real lenB, const t_blocka* at2con)
+static void assign_constraint(Lincs* li,
+ int constraint_index,
+ int a1,
+ int a2,
+ real lenA,
+ real lenB,
+ const ListOfLists<int>& at2con)
{
int con;
/* Make space in the constraint connection matrix for constraints
* connected to both end of the current constraint.
*/
- li->ncc += at2con->index[a1 + 1] - at2con->index[a1] - 1 + at2con->index[a2 + 1]
- - at2con->index[a2] - 1;
+ li->ncc += at2con[a1].ssize() - 1 + at2con[a2].ssize() - 1;
li->blnr[con + 1] = li->ncc;
/*! \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 t_blocka* at2con)
+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)
{
/* Currently this function only supports constraint groups
* in which all constraints share at least one atom
* connected constraints. We need to assign those
* to the same task.
*/
- int end;
-
- for (end = 0; end < 2; end++)
+ for (int end = 0; end < 2; end++)
{
- int a, k;
+ const int a = (end == 0 ? a1 : a2);
- a = (end == 0 ? a1 : a2);
-
- for (k = at2con->index[a]; k < at2con->index[a + 1]; k++)
+ for (const int cc : at2con[a])
{
- int cc;
-
- cc = at2con->a[k];
/* Check if constraint cc has not yet been assigned */
if (li->con_index[cc] == -1)
{
- int type;
- real lenA, lenB;
-
- type = iatom[cc * 3];
- lenA = idef.iparams[type].constr.dA;
- lenB = idef.iparams[type].constr.dB;
+ const int type = iatom[cc * 3];
+ const real lenA = idef.iparams[type].constr.dA;
+ const real lenB = idef.iparams[type].constr.dB;
if (bDynamics || lenA != 0 || lenB != 0)
{
/*! \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 t_blocka* at2con)
+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)
{
- int nca, cc[32], ca[32], k;
+ int nca, cc[32], ca[32];
int c_triangle[2] = { -1, -1 };
nca = 0;
- for (k = at2con->index[a1]; k < at2con->index[a1 + 1]; k++)
+ for (const int c : at2con[a1])
{
- int c;
-
- c = at2con->a[k];
if (c != constraint_index)
{
int aa1, aa2;
}
}
- for (k = at2con->index[a2]; k < at2con->index[a2 + 1]; k++)
+ for (const int c : at2con[a2])
{
- int c;
-
- c = at2con->a[k];
if (c != constraint_index)
{
int aa1, aa2, i;
}
//! Sets matrix indices.
-static void set_matrix_indices(Lincs* li, const Task& li_task, const t_blocka* at2con, bool bSortMatrix)
+static void set_matrix_indices(Lincs* li, const Task& li_task, const ListOfLists<int>& at2con, bool bSortMatrix)
{
for (int b = li_task.b0; b < li_task.b1; b++)
{
const int a2 = li->atoms[b].index2;
int i = li->blnr[b];
- for (int k = at2con->index[a1]; k < at2con->index[a1 + 1]; k++)
+ for (const int constraint : at2con[a1])
{
- int concon = li->con_index[at2con->a[k]];
+ const int concon = li->con_index[constraint];
if (concon != b)
{
li->blbnb[i++] = concon;
}
}
- for (int k = at2con->index[a2]; k < at2con->index[a2 + 1]; k++)
+ for (const int constraint : at2con[a2])
{
- int concon = li->con_index[at2con->a[k]];
+ const int concon = li->con_index[constraint];
if (concon != b)
{
li->blbnb[i++] = concon;
void set_lincs(const t_idef& idef, const t_mdatoms& md, bool bDynamics, const t_commrec* cr, Lincs* li)
{
int natoms;
- t_blocka at2con;
t_iatom* iatom;
li->nc_real = 0;
natoms = md.homenr;
}
- at2con = make_at2con(natoms, idef.il, idef.iparams, flexibleConstraintTreatment(bDynamics));
+ const ListOfLists<int> at2con =
+ make_at2con(natoms, idef.il, idef.iparams, flexibleConstraintTreatment(bDynamics));
const int ncon_tot = idef.il[F_CONSTR].nr / 3;
/* Skip the flexible constraints when not doing dynamics */
if (bDynamics || lenA != 0 || lenB != 0)
{
- assign_constraint(li, con, a1, a2, lenA, lenB, &at2con);
+ assign_constraint(li, con, a1, a2, lenA, lenB, at2con);
if (li->ntask > 1 && !li->bTaskDep)
{
/* We can generate independent tasks. Check if we
* need to assign connected constraints to our task.
*/
- check_assign_connected(li, iatom, idef, bDynamics, a1, a2, &at2con);
+ check_assign_connected(li, iatom, idef, bDynamics, a1, a2, at2con);
}
if (li->ntask > 1 && li->ncg_triangle > 0)
{
/* Ensure constraints in one triangle are assigned
* to the same task.
*/
- check_assign_triangle(li, iatom, idef, bDynamics, con, a1, a2, &at2con);
+ check_assign_triangle(li, iatom, idef, bDynamics, con, a1, a2, at2con);
}
}
}
li_task.tri_bits.resize(li_task.b1 - li_task.b0);
}
- set_matrix_indices(li, li_task, &at2con, bSortMatrix);
+ set_matrix_indices(li, li_task, at2con, bSortMatrix);
}
GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
}
- done_blocka(&at2con);
-
if (cr->dd == nullptr)
{
/* Since the matrix is static, we should free some memory */
struct gmx_mtop_t;
struct gmx_multisim_t;
-struct t_blocka;
struct t_commrec;
struct t_idef;
struct t_inputrec;
{
enum class ConstraintVariable : int;
-
-/* Abstract type for LINCS that is defined only in the file that uses it */
class Lincs;
+template<typename>
+class ListOfLists;
+
/*! \brief Return the data for determining constraint RMS relative deviations. */
ArrayRef<real> lincs_rmsdData(Lincs* lincsd);
real lincs_rmsd(const Lincs* lincsd);
/*! \brief Initializes and returns the lincs data struct. */
-Lincs* init_lincs(FILE* fplog,
- const gmx_mtop_t& mtop,
- int nflexcon_global,
- ArrayRef<const t_blocka> at2con,
- bool bPLINCS,
- int nIter,
- int nProjOrder);
+Lincs* init_lincs(FILE* fplog,
+ const gmx_mtop_t& mtop,
+ int nflexcon_global,
+ ArrayRef<const ListOfLists<int>> atomsToConstraintsPerMolType,
+ bool bPLINCS,
+ int nIter,
+ int nProjOrder);
/*! \brief Destructs the lincs object when it is not nullptr. */
void done_lincs(Lincs* li);
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/pbcutil/pbc.h"
-#include "gromacs/topology/block.h"
#include "gromacs/topology/idef.h"
#include "gromacs/topology/ifunc.h"
#include "gromacs/topology/topology.h"
+#include "gromacs/utility/listoflists.h"
#include "gromacs/utility/unique_cptr.h"
#include "testutils/testasserts.h"
gmx_omp_nthreads_set(emntLINCS, 1);
// Make blocka structure for faster LINCS setup
- std::vector<t_blocka> at2con_mt;
+ std::vector<ListOfLists<int>> at2con_mt;
at2con_mt.reserve(testData->mtop_.moltype.size());
for (const gmx_moltype_t& moltype : testData->mtop_.moltype)
{
&testData->nrnb_, maxwarn, &warncount_lincs);
EXPECT_TRUE(success) << "Test failed with a false return value in LINCS.";
EXPECT_EQ(warncount_lincs, 0) << "There were warnings in LINCS.";
- for (auto& moltype : at2con_mt)
- {
- sfree(moltype.index);
- sfree(moltype.a);
- }
done_lincs(lincsd);
}
#include "gromacs/math/units.h"
#include "gromacs/mdlib/constr.h"
#include "gromacs/mdtypes/inputrec.h"
-#include "gromacs/topology/block.h"
#include "gromacs/topology/idef.h"
#include "gromacs/topology/ifunc.h"
#include "gromacs/topology/topology.h"
+#include "gromacs/utility/listoflists.h"
namespace gmx
{
}
/*! \brief Returns the range of atoms constrained to atom \p a (including \p a itself) */
-static AtomIndexExtremes constraintAtomRange(int a, const t_blocka& at2con, const InteractionList& ilistConstraints)
+static AtomIndexExtremes constraintAtomRange(int a,
+ const ListOfLists<int>& at2con,
+ const InteractionList& ilistConstraints)
{
AtomIndexExtremes extremes = { a, a };
- for (int i = at2con.index[a]; i < at2con.index[a + 1]; i++)
+ for (const int constraint : at2con[a])
{
for (int j = 0; j < 2; j++)
{
- int atomJ = ilistConstraints.iatoms[at2con.a[i] * 3 + 1 + j];
+ int atomJ = ilistConstraints.iatoms[constraint * 3 + 1 + j];
extremes.minAtom = std::min(extremes.minAtom, atomJ);
extremes.maxAtom = std::max(extremes.maxAtom, atomJ);
}
}
/*! \brief Returns the size of the update group starting at \p firstAtom or 0 when criteria (see updategroups.h) are not met */
-static int detectGroup(int firstAtom,
- const gmx_moltype_t& moltype,
- const t_blocka& at2con,
- const InteractionList& ilistConstraints)
+static int detectGroup(int firstAtom,
+ const gmx_moltype_t& moltype,
+ const ListOfLists<int>& at2con,
+ const InteractionList& ilistConstraints)
{
/* We should be using moltype.atoms.atom[].ptype for checking whether
* a particle is a vsite. But the test code can't fill t_atoms,
std::vector<bool> isParticleVsite = buildIsParticleVsite(moltype);
/* A non-vsite atom without constraints is an update group by itself */
- if (!isParticleVsite[firstAtom] && at2con.index[firstAtom + 1] - at2con.index[firstAtom] == 0)
+ if (!isParticleVsite[firstAtom] && at2con[firstAtom].empty())
{
return 1;
}
}
else
{
- int numConstraints = at2con.index[a + 1] - at2con.index[a];
+ const int numConstraints = at2con[a].ssize();
if (numConstraints == 0)
{
/* We can not have unconstrained atoms in an update group */
ilistsCombined[F_CONSTR].iatoms = constraintsCombined.iatoms.data();
ilistsCombined[F_CONSTRNC].nr = 0;
/* We "include" flexible constraints, but none are present (checked above) */
- t_blocka at2con = make_at2con(moltype.atoms.nr, ilistsCombined, iparams.data(),
- FlexibleConstraintTreatment::Include);
+ const ListOfLists<int> at2con = make_at2con(moltype.atoms.nr, ilistsCombined, iparams.data(),
+ FlexibleConstraintTreatment::Include);
bool satisfiesCriteria = true;
groups.clear();
}
- done_blocka(&at2con);
-
return groups;
}
static real constraintGroupRadius(const gmx_moltype_t& moltype,
gmx::ArrayRef<const t_iparams> iparams,
const int centralAtom,
- const t_blocka& at2con,
+ const ListOfLists<int>& at2con,
const std::unordered_multimap<int, int>& angleIndices,
const real constraintLength,
const real temperature)
{
- const int numConstraints = at2con.index[centralAtom + 1] - at2con.index[centralAtom];
+ const int numConstraints = at2con[centralAtom].ssize();
GMX_RELEASE_ASSERT(numConstraints == numPartnerAtoms,
"We expect as many constraints as partner atoms here");
std::array<int, numPartnerAtoms> partnerAtoms;
for (int i = 0; i < numPartnerAtoms; i++)
{
- const int ind = at2con.a[at2con.index[centralAtom] + i] * 3;
+ const int ind = at2con[centralAtom][i] * 3;
if (ind >= moltype.ilist[F_CONSTR].size())
{
/* This is a flexible constraint, we don't optimize for that */
const InteractionList& settles = moltype.ilist[F_SETTLE];
- t_blocka at2con = make_at2con(moltype, iparams, FlexibleConstraintTreatment::Include);
+ const ListOfLists<int> at2con = make_at2con(moltype, iparams, FlexibleConstraintTreatment::Include);
const auto angleIndices = getAngleIndices(moltype);
int maxAtom = -1;
for (int a : updateGroups.block(group))
{
- int numConstraints = at2con.index[a + 1] - at2con.index[a];
+ const int numConstraints = at2con[a].ssize();
if (numConstraints > maxNumConstraints)
{
maxNumConstraints = numConstraints;
int constraintType = -1;
real maxConstraintLength = 0;
real sumConstraintLengths = 0;
- for (int i = at2con.index[maxAtom]; i < at2con.index[maxAtom + 1]; i++)
+ bool isFirstConstraint = true;
+ for (const int constraint : at2con[maxAtom])
{
- int conIndex = at2con.a[i] * (1 + NRAL(F_CONSTR));
+ int conIndex = constraint * (1 + NRAL(F_CONSTR));
int iparamsIndex;
if (conIndex < moltype.ilist[F_CONSTR].size())
{
iparamsIndex =
moltype.ilist[F_CONSTRNC].iatoms[conIndex - moltype.ilist[F_CONSTR].size()];
}
- if (i == at2con.index[maxAtom])
+ if (isFirstConstraint)
{
- constraintType = iparamsIndex;
+ constraintType = iparamsIndex;
+ isFirstConstraint = false;
}
else if (iparamsIndex != constraintType)
{
sumConstraintLengths += constraintLength;
}
- int numConstraints = at2con.index[maxAtom + 1] - at2con.index[maxAtom];
+ int numConstraints = at2con[maxAtom].ssize();
real radius;
if (numConstraints == 1)
{
maxRadius = std::max(maxRadius, dCAny);
}
- done_blocka(&at2con);
-
return maxRadius;
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff