From: Berk Hess Date: Thu, 17 May 2018 22:22:37 +0000 (+0200) Subject: Add update groups X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=cb1772dfeb0f36f073174e134f00942bfb12ee32;p=alexxy%2Fgromacs.git Add update groups Adds functionality for detecting update groups, covered by tests. Update groups are group that can be updated independently, i.e. vsite are constructed from atoms with the group and constraints only act on atoms within the group. The functionality is not used yet. Change-Id: I958983bc9010999cc7ce33b26daf587eda632c36 --- diff --git a/src/gromacs/mdlib/tests/CMakeLists.txt b/src/gromacs/mdlib/tests/CMakeLists.txt index 5bfa99cd1a..1df1d6cdb9 100644 --- a/src/gromacs/mdlib/tests/CMakeLists.txt +++ b/src/gromacs/mdlib/tests/CMakeLists.txt @@ -37,4 +37,5 @@ gmx_add_unit_test(MdlibUnitTest mdlib-test mdebin.cpp settle.cpp shake.cpp - simulationsignal.cpp) + simulationsignal.cpp + updategroups.cpp) diff --git a/src/gromacs/mdlib/tests/updategroups.cpp b/src/gromacs/mdlib/tests/updategroups.cpp new file mode 100644 index 0000000000..a52213bae7 --- /dev/null +++ b/src/gromacs/mdlib/tests/updategroups.cpp @@ -0,0 +1,280 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 2018, 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. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * http://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at http://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out http://www.gromacs.org. + */ +/*! \internal \file + * \brief + * Tests for the update groups functionality. + * + * \author berk Hess + * \ingroup module_mdlib + */ +#include "gmxpre.h" + +#include "gromacs/mdlib/updategroups.h" + +#include + +#include "gromacs/topology/topology.h" + +#include "testutils/testasserts.h" + +namespace gmx +{ + +namespace +{ + +/* TODO: Actually initialize moltype.atoms.atom when this is converted to C++ */ + +/*! \brief Returns an ethane united-atom molecule */ +gmx_moltype_t ethaneUA() +{ + gmx_moltype_t moltype = {}; + + moltype.atoms.nr = 2; + moltype.ilist[F_CONSTR].iatoms = { 0, 0, 1 }; + + return moltype; +} + +/*! \brief Returns a methane molecule */ +gmx_moltype_t methane() +{ + gmx_moltype_t moltype = {}; + + moltype.atoms.nr = 5; + moltype.ilist[F_CONSTR].iatoms = { + 0, 0, 1, + 0, 0, 2, + 0, 0, 3, + 0, 0, 4 + }; + + return moltype; +} + +/*! \brief Returns an ethane molecule */ +gmx_moltype_t ethane() +{ + gmx_moltype_t moltype = {}; + + moltype.atoms.nr = 8; + moltype.ilist[F_CONSTR].iatoms = { + 0, 0, 1, + 0, 0, 2, + 0, 0, 3, + 0, 4, 5, + 0, 4, 6, + 0, 4, 7 + }; + + return moltype; +} + +/*! \brief Returns a butane united-atom molecule */ +gmx_moltype_t butaneUA() +{ + gmx_moltype_t moltype = {}; + + moltype.atoms.nr = 4; + moltype.ilist[F_CONSTR].iatoms = { + 0, 0, 1, + 0, 1, 2, + 0, 2, 3 + }; + + return moltype; +} + +/*! \brief Returns a three-site water molecule */ +gmx_moltype_t waterThreeSite() +{ + gmx_moltype_t moltype = {}; + + moltype.atoms.nr = 3; + moltype.ilist[F_SETTLE].iatoms = { 0, 0, 1, 2}; + + return moltype; +} + +/*! \brief Returns a four-site water molecule with virtual site */ +gmx_moltype_t waterFourSite() +{ + gmx_moltype_t moltype = {}; + + moltype.atoms.nr = 4; + moltype.ilist[F_SETTLE].iatoms = { 0, 1, 2, 3 }; + moltype.ilist[F_VSITE3].iatoms = { 1, 0, 1, 2, 3}; + + return moltype; +} + +TEST(UpdateGroups, ethaneUA) +{ + gmx_mtop_t mtop; + + mtop.moltype.emplace_back(ethaneUA()); + t_iparams iparams; + iparams.constr = { 0.3, 0.3 }; + mtop.ffparams.iparams.push_back(iparams); + + auto updateGroups = gmx::makeUpdateGroups(mtop); + + ASSERT_EQ(updateGroups.size(), 1); + EXPECT_EQ(updateGroups[0].numBlocks(), 1); + + real maxRadius = computeMaxUpdateGroupRadius(mtop, updateGroups); + EXPECT_FLOAT_EQ(maxRadius, 0.3/2); +} + +TEST(UpdateGroups, methane) +{ + gmx_mtop_t mtop; + + mtop.moltype.emplace_back(methane()); + t_iparams iparams; + iparams.constr = { 0.1, 0.1 }; + mtop.ffparams.iparams.push_back(iparams); + + auto updateGroups = gmx::makeUpdateGroups(mtop); + + ASSERT_EQ(updateGroups.size(), 1); + EXPECT_EQ(updateGroups[0].numBlocks(), 1); + + real maxRadius = computeMaxUpdateGroupRadius(mtop, updateGroups); + EXPECT_FLOAT_EQ(maxRadius, 0.1); +} +TEST(UpdateGroups, ethane) +{ + gmx_mtop_t mtop; + + mtop.moltype.emplace_back(ethane()); + t_iparams iparams; + iparams.constr = { 0.1, 0.1 }; + mtop.ffparams.iparams.push_back(iparams); + + auto updateGroups = gmx::makeUpdateGroups(mtop); + + ASSERT_EQ(updateGroups.size(), 1); + EXPECT_EQ(updateGroups[0].numBlocks(), 2); + + real maxRadius = computeMaxUpdateGroupRadius(mtop, updateGroups); + EXPECT_FLOAT_EQ(maxRadius, 0.1); +} + +TEST(UpdateGroups, butaneUA) +{ + gmx_mtop_t mtop; + + mtop.moltype.emplace_back(butaneUA()); + t_iparams iparams; + iparams.constr = { 0.3, 0.3 }; + mtop.ffparams.iparams.push_back(iparams); + + auto updateGroups = gmx::makeUpdateGroups(mtop); + + EXPECT_EQ(updateGroups.size(), 0); +} + +TEST(UpdateGroups, waterThreeSite) +{ + gmx_mtop_t mtop; + + mtop.moltype.emplace_back(waterThreeSite()); + t_iparams iparams; + iparams.settle = { 0.1, 0.1633 }; + mtop.ffparams.iparams.push_back(iparams); + + auto updateGroups = gmx::makeUpdateGroups(mtop); + + ASSERT_EQ(updateGroups.size(), 1); + EXPECT_EQ(updateGroups[0].numBlocks(), 1); + + real maxRadius = computeMaxUpdateGroupRadius(mtop, updateGroups); + EXPECT_FLOAT_EQ(maxRadius, 0.083887339); +} + +// Tests update group with virtual site +TEST(UpdateGroups, waterFourSite) +{ + gmx_mtop_t mtop; + + mtop.moltype.emplace_back(waterFourSite()); + t_iparams iparams[2]; + iparams[0].settle = { 0.1, 0.1633 }; + iparams[1].vsite = { 0.128, 0.128 }; + mtop.ffparams.iparams.push_back(iparams[0]); + mtop.ffparams.iparams.push_back(iparams[1]); + + auto updateGroups = gmx::makeUpdateGroups(mtop); + + ASSERT_EQ(updateGroups.size(), 1); + EXPECT_EQ(updateGroups[0].numBlocks(), 1); +} + +TEST(UpdateGroups, fourAtomsWithSettle) +{ + gmx_mtop_t mtop; + + mtop.moltype.emplace_back(waterThreeSite()); + mtop.moltype.back().atoms.nr = 4; + + auto updateGroups = gmx::makeUpdateGroups(mtop); + + ASSERT_EQ(updateGroups.size(), 1); + EXPECT_EQ(updateGroups[0].numBlocks(), 2); +} + +TEST(UpdateGroups, twoMoltypes) +{ + gmx_mtop_t mtop; + + mtop.moltype.emplace_back(methane()); + t_iparams iparams; + iparams.constr = { 0.1, 0.1 }; + mtop.ffparams.iparams.push_back(iparams); + + mtop.moltype.emplace_back(waterThreeSite()); + // Note: iparams not accessed for SETTLE when not computing radius + + auto updateGroups = gmx::makeUpdateGroups(mtop); + + ASSERT_EQ(updateGroups.size(), 2); + EXPECT_EQ(updateGroups[0].numBlocks(), 1); + EXPECT_EQ(updateGroups[1].numBlocks(), 1); +} + +} // namespace + +} // namespace gmx diff --git a/src/gromacs/mdlib/updategroups.cpp b/src/gromacs/mdlib/updategroups.cpp new file mode 100644 index 0000000000..78309a77d3 --- /dev/null +++ b/src/gromacs/mdlib/updategroups.cpp @@ -0,0 +1,545 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 2018, 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. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * http://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at http://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out http://www.gromacs.org. + */ +/* \internal \file + * + * \brief Implements functions for generating update groups + * + * \author Berk Hess + * \ingroup module_mdlib + */ + +#include "gmxpre.h" + +#include "updategroups.h" + +#include + +#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" + +namespace gmx +{ + +/*! \brief Returns whether \p moltype contains flexible constraints */ +static bool hasFlexibleConstraints(const gmx_moltype_t &moltype, + gmx::ArrayRef iparams) +{ + for (auto &ilist : extractILists(moltype.ilist, IF_CONSTRAINT)) + { + if (ilist.functionType != F_SETTLE) + { + for (size_t i = 0; i < ilist.iatoms.size(); i += ilistStride(ilist)) + { + if (isConstraintFlexible(iparams.data(), ilist.iatoms[i])) + { + return true; + } + } + } + } + + return false; +} + +/*! \brief Returns whether moltype has incompatible vsites. + * + * 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 iparams) +{ + bool hasIncompatibleVsites = false; + + 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]]; + real coeffMin; + real coeffSum; + if (ilist.functionType == F_VSITE2) + { + coeffMin = iparam.vsite.a; + coeffSum = iparam.vsite.a; + } + else + { + coeffMin = std::min(iparam.vsite.a, iparam.vsite.b); + coeffSum = iparam.vsite.a + iparam.vsite.b; + } + if (coeffMin < 0 || coeffSum > 1) + { + hasIncompatibleVsites = true; + break; + } + } + } + else + { + hasIncompatibleVsites = true; + break; + } + } + + return hasIncompatibleVsites; +} + +/*! \brief Returns a merged list with constraints of all types */ +static InteractionList jointConstraintList(const gmx_moltype_t &moltype) +{ + InteractionList ilistCombined; + std::vector &iatoms = ilistCombined.iatoms; + + for (auto &ilist : extractILists(moltype.ilist, IF_CONSTRAINT)) + { + if (ilist.functionType == F_SETTLE) + { + for (size_t i = 0; i < ilist.iatoms.size(); i += ilistStride(ilist)) + { + iatoms.push_back(-1); + iatoms.push_back(ilist.iatoms[i + 1]); + iatoms.push_back(ilist.iatoms[i + 2]); + iatoms.push_back(-1); + iatoms.push_back(ilist.iatoms[i + 1]); + iatoms.push_back(ilist.iatoms[i + 3]); + iatoms.push_back(-1); + iatoms.push_back(ilist.iatoms[i + 2]); + iatoms.push_back(ilist.iatoms[i + 3]); + } + } + else + { + 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()); + } + } + + return ilistCombined; +} + +/*! \brief Struct for returning an atom range */ +struct AtomIndexExtremes +{ + int minAtom; //!< The minimum atom index + int maxAtom; //!< The maximum atom index +}; + +/*! \brief Returns the range of constructing atom for vsite with atom index \p a */ +static AtomIndexExtremes +vsiteConstructRange(int a, + const gmx_moltype_t &moltype) +{ + AtomIndexExtremes extremes = { -1, -1 }; + + for (auto &ilist : extractILists(moltype.ilist, IF_VSITE)) + { + for (size_t i = 0; i < ilist.iatoms.size(); i += ilistStride(ilist)) + { + if (ilist.iatoms[i + 1] == a) + { + extremes.minAtom = ilist.iatoms[i + 2]; + extremes.maxAtom = ilist.iatoms[i + 2]; + for (size_t j = i + 3; j < i + ilistStride(ilist); j++) + { + extremes.minAtom = std::min(extremes.minAtom, ilist.iatoms[j]); + extremes.maxAtom = std::max(extremes.maxAtom, ilist.iatoms[j]); + } + return extremes; + } + } + } + + GMX_RELEASE_ASSERT(false, "If a is a vsite, we should have found constructing atoms"); + + return extremes; +} + +/*! \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) +{ + AtomIndexExtremes extremes = { a, a }; + + for (int i = at2con.index[a]; i < at2con.index[a + 1]; i++) + { + for (int j = 0; j < 2; 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); + } + } + + return extremes; +} + +/*! \brief Returns a list that tells whether atoms in \p moltype are vsites */ +static std::vector buildIsParticleVsite(const gmx_moltype_t &moltype) +{ + std::vector isVsite(moltype.atoms.nr); + + for (auto &ilist : extractILists(moltype.ilist, IF_VSITE)) + { + for (size_t i = 0; i < ilist.iatoms.size(); i += ilistStride(ilist)) + { + int vsiteAtom = ilist.iatoms[i + 1]; + isVsite[vsiteAtom] = true; + } + } + + return isVsite; +} + +/*! \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) +{ + /* 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, + * because it uses C pointers which get double freed. + */ + std::vector 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) + { + return 1; + } + + /* A (potential) update group starts at firstAtom and should consist + * of two or more atoms and possibly vsites. At least one atom should + * have constraints with all other atoms and vsites should have all + * constructing atoms inside the group. Here we increase lastAtom until + * the criteria are fulfilled or exit when criteria cannot be met. + */ + int numAtomsWithConstraints = 0; + int maxConstraintsPerAtom = 0; + int lastAtom = firstAtom; + int a = firstAtom; + while (a <= lastAtom) + { + if (isParticleVsite[a]) + { + AtomIndexExtremes extremes = vsiteConstructRange(a, moltype); + if (extremes.minAtom < firstAtom) + { + /* A constructing atom is outside the group, + * we can not use update groups. + */ + return 0; + } + lastAtom = std::max(lastAtom, extremes.maxAtom); + } + else + { + int numConstraints = at2con.index[a + 1] - at2con.index[a]; + if (numConstraints == 0) + { + /* We can not have unconstrained atoms in an update group */ + return 0; + } + /* This atom has at least one constraint. + * Check whether all constraints are within the group + * and whether we need to extend the group. + */ + numAtomsWithConstraints += 1; + maxConstraintsPerAtom = std::max(maxConstraintsPerAtom, numConstraints); + + AtomIndexExtremes extremes = + constraintAtomRange(a, at2con, ilistConstraints); + if (extremes.minAtom < firstAtom) + { + /* Constraint to atom outside the "group" */ + return 0; + } + lastAtom = std::max(lastAtom, extremes.maxAtom); + } + + a++; + } + + /* 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]) + { + AtomIndexExtremes extremes = vsiteConstructRange(lastAtom + 1, moltype); + if (extremes.minAtom < firstAtom) + { + /* Constructing atom precedes the group */ + return 0; + } + else if (extremes.maxAtom <= lastAtom) + { + /* All constructing atoms are in the group, add the vsite to the group */ + lastAtom++; + } + else if (extremes.minAtom <= lastAtom) + { + /* Some, but not all constructing atoms are in the group */ + return 0; + } + } + + GMX_RELEASE_ASSERT(maxConstraintsPerAtom < numAtomsWithConstraints, + "We have checked that atoms are only constrained to atoms within the group," + "so each atom should have fewer constraints than the number of atoms"); + /* Check that at least one atom is constrained to all others */ + if (maxConstraintsPerAtom != numAtomsWithConstraints - 1) + { + return 0; + } + + return lastAtom - firstAtom + 1; +} + +/*! \brief Returns a list of update groups for \p moltype */ +static RangePartitioning +makeUpdateGroups(const gmx_moltype_t &moltype, + gmx::ArrayRef iparams) +{ + RangePartitioning groups; + + /* Update groups are not compatible with flexible constraints. + * Without dynamics the flexible constraints are ignored, + * 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)) + { + 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; + /* We "include" flexible constraints, but none are present (checked above) */ + t_blocka at2con = make_at2con(moltype.atoms.nr, + ilistsCombined, iparams.data(), + FlexibleConstraintTreatment::Include); + + bool satisfiesCriteria = true; + + int firstAtom = 0; + while (satisfiesCriteria && firstAtom < moltype.atoms.nr) + { + int numAtomsInGroup = + detectGroup(firstAtom, moltype, at2con, constraintsCombined); + + if (numAtomsInGroup == 0) + { + satisfiesCriteria = false; + } + else + { + groups.appendBlock(numAtomsInGroup); + } + firstAtom += numAtomsInGroup; + } + + if (!satisfiesCriteria) + { + /* Make groups empty, to signal not satisfying the criteria */ + groups.clear(); + } + + done_blocka(&at2con); + + return groups; +} + +std::vector makeUpdateGroups(const gmx_mtop_t &mtop) +{ + std::vector updateGroups; + + bool systemSatisfiesCriteria = true; + for (const gmx_moltype_t &moltype : mtop.moltype) + { + updateGroups.push_back(makeUpdateGroups(moltype, mtop.ffparams.iparams)); + + if (updateGroups.back().numBlocks() == 0) + { + systemSatisfiesCriteria = false; + } + } + + if (!systemSatisfiesCriteria) + { + updateGroups.clear(); + } + + return updateGroups; +} + +/*! \brief Returns the maximum update group radius for \p moltype */ +static real +computeMaxUpdateGroupRadius(const gmx_moltype_t &moltype, + gmx::ArrayRef iparams, + const RangePartitioning &updateGroups) +{ + GMX_RELEASE_ASSERT(!hasFlexibleConstraints(moltype, iparams), + "Flexible constraints are not supported here"); + + const InteractionList &settles = moltype.ilist[F_SETTLE]; + + t_blocka at2con = + make_at2con(moltype, iparams, FlexibleConstraintTreatment::Include); + + const int stride = 1 + NRAL(F_CONSTR); + const real half = 0.5; + + real maxRadius = 0; + for (int group = 0; group < updateGroups.numBlocks(); group++) + { + if (updateGroups.block(group).size() == 1) + { + /* Single atom group, radius is zero */ + continue; + } + + /* Find the atom maxAtom with the maximum number of constraints */ + int maxNumConstraints = 0; + int maxAtom = -1; + for (int a : updateGroups.block(group)) + { + int numConstraints = at2con.index[a + 1] - at2con.index[a]; + if (numConstraints > maxNumConstraints) + { + maxNumConstraints = numConstraints; + maxAtom = a; + } + } + GMX_ASSERT(maxAtom >= 0 || settles.size() > 0, + "We should have at least two atoms in the group with constraints"); + if (maxAtom < 0) + { + continue; + } + + for (int i = at2con.index[maxAtom]; i < at2con.index[maxAtom + 1]; i++) + { + int conIndex = at2con.a[i]*stride; + int iparamsIndex; + if (conIndex < moltype.ilist[F_CONSTR].size()) + { + iparamsIndex = moltype.ilist[F_CONSTR].iatoms[conIndex]; + } + else + { + iparamsIndex = moltype.ilist[F_CONSTRNC].iatoms[conIndex - moltype.ilist[F_CONSTR].size()]; + } + /* Here we take the maximum constraint length of the A and B + * 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); + if (at2con.index[maxAtom + 1] == at2con.index[maxAtom] + 1) + { + /* Single constraint: use the distance from the midpoint */ + maxRadius = std::max(maxRadius, half*constraintLength); + } + else + { + /* Multiple constraints: use the distance from the central + * atom. We can do better than this if we would assume + * that the angles between bonds do not stretch a lot. + */ + maxRadius = std::max(maxRadius, constraintLength); + } + } + } + + 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; + /* 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 dCAny = std::sqrt(std::max(dCO2, dCH2)); + maxRadius = std::max(maxRadius, dCAny); + } + + done_blocka(&at2con); + + return maxRadius; +} + +real computeMaxUpdateGroupRadius(const gmx_mtop_t &mtop, + gmx::ArrayRef updateGroups) +{ + if (updateGroups.empty()) + { + return 0; + } + + GMX_RELEASE_ASSERT(static_cast(updateGroups.size()) == mtop.moltype.size(), + "We need one update group entry per moleculetype"); + + real maxRadius = 0; + + for (size_t moltype = 0; moltype < mtop.moltype.size(); moltype++) + { + maxRadius + = std::max(maxRadius, + computeMaxUpdateGroupRadius(mtop.moltype[moltype], + mtop.ffparams.iparams, + updateGroups[moltype])); + } + + return maxRadius; +} + +} // namespace gmx diff --git a/src/gromacs/mdlib/updategroups.h b/src/gromacs/mdlib/updategroups.h new file mode 100644 index 0000000000..6ab9c35dd0 --- /dev/null +++ b/src/gromacs/mdlib/updategroups.h @@ -0,0 +1,88 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 2018, 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. + * + * GROMACS is free software; you can redistribute it and/or + * modify it under the terms of the GNU Lesser General Public License + * as published by the Free Software Foundation; either version 2.1 + * of the License, or (at your option) any later version. + * + * GROMACS is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + * Lesser General Public License for more details. + * + * You should have received a copy of the GNU Lesser General Public + * License along with GROMACS; if not, see + * http://www.gnu.org/licenses, or write to the Free Software Foundation, + * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA. + * + * If you want to redistribute modifications to GROMACS, please + * consider that scientific software is very special. Version + * control is crucial - bugs must be traceable. We will be happy to + * consider code for inclusion in the official distribution, but + * derived work must not be called official GROMACS. Details are found + * in the README & COPYING files - if they are missing, get the + * official version at http://www.gromacs.org. + * + * To help us fund GROMACS development, we humbly ask that you cite + * the research papers on the package. Check out http://www.gromacs.org. + */ +/*! \libinternal \file + * + * \brief Declares the functions for generating update groups + * + * \author Berk Hess + * \ingroup module_mdlib + * \inlibraryapi + */ +#ifndef GMX_MDLIB_UPDATEGROUPS +#define GMX_MDLIB_UPDATEGROUPS + +#include + +#include "gromacs/utility/arrayref.h" +#include "gromacs/utility/real.h" + +struct gmx_mtop_t; +struct t_inputrec; + +namespace gmx +{ +class RangePartitioning; + +/*! \brief Returns a vector with update groups for each moleculetype in \p mtop + * or an empty vector when the criteria (see below) are not satisfied. + * + * An empty vector is returned when at least one moleculetype does not obey + * the restrictions of update groups, e.g. more than two constraints in a row. + * + * Currently valid update groups are: + * a single atom which is not a virtual site and does not have constraints; + * or a group of atoms where all virtual sites are constructed from atoms + * within the group and at least one non-vsite atom is constrained to + * all other non-vsite atoms. + * To have update groups, all virtual sites should be linear 2 or 3 atom + * constructions with coefficients >= 0 and sum of coefficients <= 1. + * + * \param[in] mtop The system topology + */ +std::vector makeUpdateGroups(const gmx_mtop_t &mtop); + +/*! \brief Returns the maximum update group radius + * + * \note When \p updateGroups is empty, 0 is returned. + * + * \param[in] mtop The system topology + * \param[in] updateGroups List of update group, size should match the number of moltypes in \p mtop or be 0 + */ +real computeMaxUpdateGroupRadius(const gmx_mtop_t &mtop, + gmx::ArrayRef updateGroups); + +} // namespace gmx + +#endif // GMX_MDLIB_UPDATEGROUPS diff --git a/src/gromacs/topology/idef.h b/src/gromacs/topology/idef.h index 7c4f536ce2..62a9fa6f0d 100644 --- a/src/gromacs/topology/idef.h +++ b/src/gromacs/topology/idef.h @@ -271,6 +271,15 @@ extractILists(const InteractionLists &ilists, return handles; } +/*! \brief Returns the stride for the iatoms array in \p ilistHandle + * + * \param[in] ilistHandle The ilist to return the stride for + */ +static inline int ilistStride(const InteractionListHandle &ilistHandle) +{ + return 1 + NRAL(ilistHandle.functionType); +} + struct gmx_cmapdata_t { std::vector cmap; /* Has length 4*grid_spacing*grid_spacing, */