{
/*! \brief Returns whether \p moltype contains flexible constraints */
-static bool hasFlexibleConstraints(const gmx_moltype_t &moltype,
- gmx::ArrayRef<const t_iparams> iparams)
+static bool hasFlexibleConstraints(const gmx_moltype_t& moltype, gmx::ArrayRef<const t_iparams> iparams)
{
- for (auto &ilist : extractILists(moltype.ilist, IF_CONSTRAINT))
+ for (auto& ilist : extractILists(moltype.ilist, IF_CONSTRAINT))
{
if (ilist.functionType != F_SETTLE)
{
* For simplicity the only compatible vsites are linear 2 or 3 atom sites
* that are constructed in between the 2 or 3 contructing atoms,
*/
-static bool hasIncompatibleVsites(const gmx_moltype_t &moltype,
- gmx::ArrayRef<const t_iparams> iparams)
+static bool hasIncompatibleVsites(const gmx_moltype_t& moltype, gmx::ArrayRef<const t_iparams> iparams)
{
bool hasIncompatibleVsites = false;
- for (auto &ilist : extractILists(moltype.ilist, IF_VSITE))
+ for (auto& ilist : extractILists(moltype.ilist, IF_VSITE))
{
if (ilist.functionType == F_VSITE2 || ilist.functionType == F_VSITE3)
{
for (size_t i = 0; i < ilist.iatoms.size(); i += ilistStride(ilist))
{
- const t_iparams &iparam = iparams[ilist.iatoms[i]];
+ const t_iparams& iparam = iparams[ilist.iatoms[i]];
real coeffMin;
real coeffSum;
if (ilist.functionType == F_VSITE2)
}
/*! \brief Returns a merged list with constraints of all types */
-static InteractionList jointConstraintList(const gmx_moltype_t &moltype)
+static InteractionList jointConstraintList(const gmx_moltype_t& moltype)
{
InteractionList ilistCombined;
- std::vector<int> &iatoms = ilistCombined.iatoms;
+ std::vector<int>& iatoms = ilistCombined.iatoms;
- for (auto &ilist : extractILists(moltype.ilist, IF_CONSTRAINT))
+ for (auto& ilist : extractILists(moltype.ilist, IF_CONSTRAINT))
{
if (ilist.functionType == F_SETTLE)
{
}
else
{
- GMX_RELEASE_ASSERT(NRAL(ilist.functionType) == 2, "Can only handle two-atom non-SETTLE constraints");
+ GMX_RELEASE_ASSERT(NRAL(ilist.functionType) == 2,
+ "Can only handle two-atom non-SETTLE constraints");
- iatoms.insert(iatoms.end(),
- ilist.iatoms.begin(), ilist.iatoms.end());
+ iatoms.insert(iatoms.end(), ilist.iatoms.begin(), ilist.iatoms.end());
}
}
};
/*! \brief Returns the range of constructing atom for vsite with atom index \p a */
-static AtomIndexExtremes
-vsiteConstructRange(int a,
- const gmx_moltype_t &moltype)
+static AtomIndexExtremes vsiteConstructRange(int a, const gmx_moltype_t& moltype)
{
AtomIndexExtremes extremes = { -1, -1 };
- for (auto &ilist : extractILists(moltype.ilist, IF_VSITE))
+ for (auto& ilist : extractILists(moltype.ilist, IF_VSITE))
{
for (size_t i = 0; i < ilist.iatoms.size(); i += ilistStride(ilist))
{
}
/*! \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 t_blocka& at2con, const InteractionList& ilistConstraints)
{
AtomIndexExtremes extremes = { a, a };
{
for (int j = 0; j < 2; j++)
{
- int atomJ = ilistConstraints.iatoms[at2con.a[i]*3 + 1 + j];
+ int atomJ = ilistConstraints.iatoms[at2con.a[i] * 3 + 1 + j];
extremes.minAtom = std::min(extremes.minAtom, atomJ);
extremes.maxAtom = std::max(extremes.maxAtom, atomJ);
}
}
/*! \brief Returns a list that tells whether atoms in \p moltype are vsites */
-static std::vector<bool> buildIsParticleVsite(const gmx_moltype_t &moltype)
+static std::vector<bool> buildIsParticleVsite(const gmx_moltype_t& moltype)
{
std::vector<bool> isVsite(moltype.atoms.nr);
- for (auto &ilist : extractILists(moltype.ilist, IF_VSITE))
+ for (auto& ilist : extractILists(moltype.ilist, IF_VSITE))
{
for (size_t i = 0; i < ilist.iatoms.size(); i += ilistStride(ilist))
{
/*! \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)
+ const gmx_moltype_t& moltype,
+ const t_blocka& 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.index[firstAtom + 1] - at2con.index[firstAtom] == 0)
{
return 1;
}
* and whether we need to extend the group.
*/
numAtomsWithConstraints += 1;
- maxConstraintsPerAtom = std::max(maxConstraintsPerAtom, numConstraints);
+ maxConstraintsPerAtom = std::max(maxConstraintsPerAtom, numConstraints);
- AtomIndexExtremes extremes =
- constraintAtomRange(a, at2con, ilistConstraints);
+ AtomIndexExtremes extremes = constraintAtomRange(a, at2con, ilistConstraints);
if (extremes.minAtom < firstAtom)
{
/* Constraint to atom outside the "group" */
/* lastAtom might be followed by a vsite that is constructed from atoms
* with index <= lastAtom. Handle this case.
*/
- if (lastAtom + 1 < moltype.atoms.nr &&
- isParticleVsite[lastAtom + 1])
+ if (lastAtom + 1 < moltype.atoms.nr && isParticleVsite[lastAtom + 1])
{
AtomIndexExtremes extremes = vsiteConstructRange(lastAtom + 1, moltype);
if (extremes.minAtom < firstAtom)
}
/*! \brief Returns a list of update groups for \p moltype */
-static RangePartitioning
-makeUpdateGroups(const gmx_moltype_t &moltype,
- gmx::ArrayRef<const t_iparams> iparams)
+static RangePartitioning makeUpdateGroups(const gmx_moltype_t& moltype, gmx::ArrayRef<const t_iparams> iparams)
{
RangePartitioning groups;
* but since performance for EM/NM is less critical, we do not
* use update groups to keep the code here simpler.
*/
- if (hasFlexibleConstraints(moltype, iparams) ||
- hasIncompatibleVsites(moltype, iparams))
+ if (hasFlexibleConstraints(moltype, iparams) || hasIncompatibleVsites(moltype, iparams))
{
return groups;
}
/* 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;
+ ilistsCombined[F_CONSTR].nr = constraintsCombined.iatoms.size();
+ 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);
+ t_blocka at2con = make_at2con(moltype.atoms.nr, ilistsCombined, iparams.data(),
+ FlexibleConstraintTreatment::Include);
bool satisfiesCriteria = true;
- int firstAtom = 0;
+ int firstAtom = 0;
while (satisfiesCriteria && firstAtom < moltype.atoms.nr)
{
- int numAtomsInGroup =
- detectGroup(firstAtom, moltype, at2con, constraintsCombined);
+ int numAtomsInGroup = detectGroup(firstAtom, moltype, at2con, constraintsCombined);
if (numAtomsInGroup == 0)
{
return groups;
}
-std::vector<RangePartitioning> makeUpdateGroups(const gmx_mtop_t &mtop)
+std::vector<RangePartitioning> makeUpdateGroups(const gmx_mtop_t& mtop)
{
std::vector<RangePartitioning> updateGroups;
- bool systemSatisfiesCriteria = true;
- for (const gmx_moltype_t &moltype : mtop.moltype)
+ bool systemSatisfiesCriteria = true;
+ for (const gmx_moltype_t& moltype : mtop.moltype)
{
updateGroups.push_back(makeUpdateGroups(moltype, mtop.ffparams.iparams));
}
/*! \brief Returns a map of angles ilist.iatoms indices with the middle atom as key */
-static std::unordered_multimap<int, int>
-getAngleIndices(const gmx_moltype_t &moltype)
+static std::unordered_multimap<int, int> getAngleIndices(const gmx_moltype_t& moltype)
{
- const InteractionList &angles = moltype.ilist[F_ANGLES];
+ const InteractionList& angles = moltype.ilist[F_ANGLES];
std::unordered_multimap<int, int> indices(angles.size());
* by one angle potential, when a potential is perturbed or when an angle
* could reach more than 180 degrees.
*/
-template <int numPartnerAtoms> static real
-constraintGroupRadius(const gmx_moltype_t &moltype,
- gmx::ArrayRef<const t_iparams> iparams,
- const int centralAtom,
- const t_blocka &at2con,
- const std::unordered_multimap<int, int> &angleIndices,
- const real constraintLength,
- const real temperature)
+template<int numPartnerAtoms>
+static real constraintGroupRadius(const gmx_moltype_t& moltype,
+ gmx::ArrayRef<const t_iparams> iparams,
+ const int centralAtom,
+ const t_blocka& at2con,
+ const std::unordered_multimap<int, int>& angleIndices,
+ const real constraintLength,
+ const real temperature)
{
const int numConstraints = at2con.index[centralAtom + 1] - at2con.index[centralAtom];
- GMX_RELEASE_ASSERT(numConstraints == numPartnerAtoms, "We expect as many constraints as partner atoms here");
+ 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.a[at2con.index[centralAtom] + i] * 3;
if (ind >= moltype.ilist[F_CONSTR].size())
{
/* This is a flexible constraint, we don't optimize for that */
partnerAtoms[i] = (a1 == centralAtom ? a2 : a1);
}
- const InteractionList &angles = moltype.ilist[F_ANGLES];
- auto range = angleIndices.equal_range(centralAtom);
- int angleType = -1;
- std::array<int, numPartnerAtoms> numAngles = { 0 };
- bool areSameType = true;
+ const InteractionList& angles = moltype.ilist[F_ANGLES];
+ auto range = angleIndices.equal_range(centralAtom);
+ int angleType = -1;
+ std::array<int, numPartnerAtoms> numAngles = { 0 };
+ bool areSameType = true;
for (auto it = range.first; it != range.second; ++it)
{
/* Check if the outer atoms in the angle are both partner atoms */
int numAtomsFound = 0;
for (int ind = it->second + 1; ind < it->second + 4; ind += 2)
{
- for (const int &partnerAtom : partnerAtoms)
+ for (const int& partnerAtom : partnerAtoms)
{
if (angles.iatoms[ind] == partnerAtom)
{
}
/* We don't bother optimizing the perturbed angle case */
- const t_iparams &angleParams = iparams[angleType];
- if (criteriaSatisfied &&
- angleParams.harmonic.rB == angleParams.harmonic.rA &&
- angleParams.harmonic.krB == angleParams.harmonic.krA)
+ const t_iparams& angleParams = iparams[angleType];
+ if (criteriaSatisfied && angleParams.harmonic.rB == angleParams.harmonic.rA
+ && angleParams.harmonic.krB == angleParams.harmonic.krA)
{
/* Set number of stddevs such that change of exceeding < 10^-9 */
constexpr real c_numSigma = 6.0;
/* Compute the maximally stretched angle */
- const real eqAngle = angleParams.harmonic.rA*DEG2RAD;
- const real fc = angleParams.harmonic.krA;
- const real maxAngle = eqAngle + c_numSigma*BOLTZ*temperature/((numPartnerAtoms - 1)*fc);
+ const real eqAngle = angleParams.harmonic.rA * DEG2RAD;
+ const real fc = angleParams.harmonic.krA;
+ const real maxAngle = eqAngle + c_numSigma * BOLTZ * temperature / ((numPartnerAtoms - 1) * fc);
if (maxAngle >= M_PI)
{
return -1;
* Return the distance from the COG to the farthest two corners,
* i.e. the partner atoms.
*/
- real distMidPartner = std::sin(0.5*maxAngle)*constraintLength;
- real distCentralMid = std::cos(0.5*maxAngle)*constraintLength;
- real distCogMid = distCentralMid*numPartnerAtoms/(numPartnerAtoms + 1);
+ real distMidPartner = std::sin(0.5 * maxAngle) * constraintLength;
+ real distCentralMid = std::cos(0.5 * maxAngle) * constraintLength;
+ real distCogMid = distCentralMid * numPartnerAtoms / (numPartnerAtoms + 1);
real distCogPartner = std::sqrt(gmx::square(distMidPartner) + gmx::square(distCogMid));
return distCogPartner;
* Then we compute the distance of the central atom to the plane.
* The distance of the COG to the ourlier is returned.
*/
- real halfDistEqPartner = std::sin(0.5*eqAngle)*constraintLength;
- real distPartnerOutlier = 2*std::sin(0.5*maxAngle)*constraintLength;
- real distMidOutlier = std::sqrt(gmx::square(distPartnerOutlier) - gmx::square(halfDistEqPartner));
- real distMidCenterInPlane = 0.5*(distMidOutlier - gmx::square(halfDistEqPartner)/distMidOutlier);
- real distCentralToPlane = std::sqrt(gmx::square(constraintLength) - gmx::square(halfDistEqPartner) - gmx::square(distMidCenterInPlane));
- real distCogOutlierH = distCentralToPlane/(numPartnerAtoms + 1);
- real distCogOutlierP = distMidOutlier - (distMidOutlier + distMidCenterInPlane)/(numPartnerAtoms + 1);
- real distCogOutlier = std::sqrt(gmx::square(distCogOutlierH) + gmx::square(distCogOutlierP));
+ real halfDistEqPartner = std::sin(0.5 * eqAngle) * constraintLength;
+ real distPartnerOutlier = 2 * std::sin(0.5 * maxAngle) * constraintLength;
+ real distMidOutlier =
+ std::sqrt(gmx::square(distPartnerOutlier) - gmx::square(halfDistEqPartner));
+ real distMidCenterInPlane =
+ 0.5 * (distMidOutlier - gmx::square(halfDistEqPartner) / distMidOutlier);
+ real distCentralToPlane =
+ std::sqrt(gmx::square(constraintLength) - gmx::square(halfDistEqPartner)
+ - gmx::square(distMidCenterInPlane));
+ real distCogOutlierH = distCentralToPlane / (numPartnerAtoms + 1);
+ real distCogOutlierP =
+ distMidOutlier - (distMidOutlier + distMidCenterInPlane) / (numPartnerAtoms + 1);
+ real distCogOutlier = std::sqrt(gmx::square(distCogOutlierH) + gmx::square(distCogOutlierP));
return distCogOutlier;
}
}
/*! \brief Returns the maximum update group radius for \p moltype */
-static real
-computeMaxUpdateGroupRadius(const gmx_moltype_t &moltype,
- gmx::ArrayRef<const t_iparams> iparams,
- const RangePartitioning &updateGroups,
- real temperature)
+static real computeMaxUpdateGroupRadius(const gmx_moltype_t& moltype,
+ gmx::ArrayRef<const t_iparams> iparams,
+ const RangePartitioning& updateGroups,
+ real temperature)
{
GMX_RELEASE_ASSERT(!hasFlexibleConstraints(moltype, iparams),
"Flexible constraints are not supported here");
- const InteractionList &settles = moltype.ilist[F_SETTLE];
+ const InteractionList& settles = moltype.ilist[F_SETTLE];
- t_blocka at2con =
- make_at2con(moltype, iparams, FlexibleConstraintTreatment::Include);
+ t_blocka at2con = make_at2con(moltype, iparams, FlexibleConstraintTreatment::Include);
const auto angleIndices = getAngleIndices(moltype);
- real maxRadius = 0;
+ real maxRadius = 0;
for (int group = 0; group < updateGroups.numBlocks(); group++)
{
if (updateGroups.block(group).size() == 1)
real sumConstraintLengths = 0;
for (int i = at2con.index[maxAtom]; i < at2con.index[maxAtom + 1]; i++)
{
- int conIndex = at2con.a[i]*(1 + NRAL(F_CONSTR));
+ int conIndex = at2con.a[i] * (1 + NRAL(F_CONSTR));
int iparamsIndex;
if (conIndex < moltype.ilist[F_CONSTR].size())
{
}
else
{
- iparamsIndex = moltype.ilist[F_CONSTRNC].iatoms[conIndex - moltype.ilist[F_CONSTR].size()];
+ iparamsIndex =
+ moltype.ilist[F_CONSTRNC].iatoms[conIndex - moltype.ilist[F_CONSTR].size()];
}
if (i == at2con.index[maxAtom])
{
* topology, which assumes lambda is between 0 and 1 for
* free-energy runs.
*/
- real constraintLength = std::max(iparams[iparamsIndex].constr.dA,
- iparams[iparamsIndex].constr.dB);
- maxConstraintLength = std::max(maxConstraintLength,
- constraintLength);
- sumConstraintLengths += constraintLength;
+ real constraintLength =
+ std::max(iparams[iparamsIndex].constr.dA, iparams[iparamsIndex].constr.dB);
+ maxConstraintLength = std::max(maxConstraintLength, constraintLength);
+ sumConstraintLengths += constraintLength;
}
int numConstraints = at2con.index[maxAtom + 1] - at2con.index[maxAtom];
if (numConstraints == 1)
{
/* Single constraint: the radius is the distance from the midpoint */
- radius = 0.5_real*maxConstraintLength;
+ radius = 0.5_real * maxConstraintLength;
}
else
{
*/
if (numConstraints == 2 && allTypesAreEqual && temperature > 0)
{
- radius = constraintGroupRadius<2>(moltype, iparams,
- maxAtom, at2con,
- angleIndices,
- maxConstraintLength,
- temperature);
+ radius = constraintGroupRadius<2>(moltype, iparams, maxAtom, at2con, angleIndices,
+ maxConstraintLength, temperature);
}
/* With 3 constraints the maximum possible radius is 1.4 times
* the constraint length, so it is worth computing a smaller
*/
if (numConstraints == 3 && allTypesAreEqual && temperature >= 0)
{
- radius = constraintGroupRadius<3>(moltype, iparams,
- maxAtom, at2con,
- angleIndices,
- maxConstraintLength,
- temperature);
+ radius = constraintGroupRadius<3>(moltype, iparams, maxAtom, at2con, angleIndices,
+ maxConstraintLength, temperature);
if (temperature == 0 && radius >= 0)
{
/* Add a 10% margin for deviation at 0 K */
/* Worst case: atom with the longest constraint on one side
* of the center, all others on the opposite side
*/
- radius = maxConstraintLength + (sumConstraintLengths - 2*maxConstraintLength)/(numConstraints + 1);
+ radius = maxConstraintLength
+ + (sumConstraintLengths - 2 * maxConstraintLength) / (numConstraints + 1);
}
}
maxRadius = std::max(maxRadius, radius);
for (int i = 0; i < settles.size(); i += 1 + NRAL(F_SETTLE))
{
- const real dOH = iparams[settles.iatoms[i]].settle.doh;
- const real dHH = iparams[settles.iatoms[i]].settle.dhh;
+ const real dOH = iparams[settles.iatoms[i]].settle.doh;
+ const real dHH = iparams[settles.iatoms[i]].settle.dhh;
/* Compute distance^2 from center of geometry to O and H */
- const real dCO2 = (4*dOH*dOH - dHH*dHH)/9;
- const real dCH2 = (dOH*dOH + 2*dHH*dHH)/9;
+ const real dCO2 = (4 * dOH * dOH - dHH * dHH) / 9;
+ const real dCH2 = (dOH * dOH + 2 * dHH * dHH) / 9;
const real dCAny = std::sqrt(std::max(dCO2, dCH2));
maxRadius = std::max(maxRadius, dCAny);
}
return maxRadius;
}
-real
-computeMaxUpdateGroupRadius(const gmx_mtop_t &mtop,
- gmx::ArrayRef<const RangePartitioning> updateGroups,
- real temperature)
+real computeMaxUpdateGroupRadius(const gmx_mtop_t& mtop,
+ gmx::ArrayRef<const RangePartitioning> updateGroups,
+ real temperature)
{
if (updateGroups.empty())
{
for (size_t moltype = 0; moltype < mtop.moltype.size(); moltype++)
{
- maxRadius
- = std::max(maxRadius,
- computeMaxUpdateGroupRadius(mtop.moltype[moltype],
- mtop.ffparams.iparams,
- updateGroups[moltype],
- temperature));
+ maxRadius = std::max(
+ maxRadius, computeMaxUpdateGroupRadius(mtop.moltype[moltype], mtop.ffparams.iparams,
+ updateGroups[moltype], temperature));
}
return maxRadius;
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff