Make index groups available for mdrun modules during preprocessing.
Notify mdrun modules with a DefaultIndexGroupsAndNames class that
translates index group names into atom indices.
refs #2282
Change-Id: Idc1bc3a5840e7b7ec3c82904416887b02417c92c
fprintf(stderr, "initialising group options...\n");
}
do_index(mdparin, ftp2fn_null(efNDX, NFILE, fnm),
- &sys, bVerbose, ir,
- wi);
+ &sys, bVerbose, mdModules.notifier(), ir, wi);
if (ir->cutoff_scheme == ecutsVERLET && ir->verletbuf_tol > 0)
{
#include "gromacs/options/options.h"
#include "gromacs/options/treesupport.h"
#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/selection/indexutil.h"
#include "gromacs/topology/block.h"
#include "gromacs/topology/ifunc.h"
#include "gromacs/topology/index.h"
void do_index(const char* mdparin, const char *ndx,
gmx_mtop_t *mtop,
bool bVerbose,
+ const gmx::MDModules::notifier_type ¬ifier,
t_inputrec *ir,
warninp_t wi)
{
make_IMD_group(ir->imd, is->imd_grp, defaultIndexGroups, gnames);
}
+ gmx::IndexGroupsAndNames defaultIndexGroupsAndNames(
+ *defaultIndexGroups, gmx::arrayRefFromArray(gnames, defaultIndexGroups->nr));
+ notifier.notify(defaultIndexGroupsAndNames);
+
auto accelerations = gmx::splitString(is->acc);
auto accelerationGroupNames = gmx::splitString(is->accgrps);
if (accelerationGroupNames.size() * DIM != accelerations.size())
#include "gromacs/fileio/readinp.h"
#include "gromacs/math/vectypes.h"
+#include "gromacs/mdrun/mdmodules.h"
#include "gromacs/utility/real.h"
namespace gmx
* function is called. Also prints the input file back to mdparout.
*/
-void do_index(const char* mdparin,
- const char *ndx,
- gmx_mtop_t *mtop,
- bool bVerbose,
- t_inputrec *ir,
- warninp_t wi);
+void do_index(const char * mdparin,
+ const char *ndx,
+ gmx_mtop_t *mtop,
+ bool bVerbose,
+ const gmx::MDModules::notifier_type ¬ifier,
+ t_inputrec *ir,
+ warninp_t wi);
/* Read the index file and assign grp numbers to atoms.
*/
class KeyValueTreeBuilder;
class IMDModule;
class LocalAtomSetManager;
+class IndexGroupsAndNames;
/*! \libinternal \brief
* \brief Signals that the communication record is set up and provides this record.
//! Register callback function types for MDModules
using notifier_type = registerMdModuleNotification<
CommunicationIsSetup,
+ IndexGroupsAndNames,
KeyValueTreeBuilder*,
const KeyValueTreeObject &,
LocalAtomSetManager *
#include <cstring>
#include <algorithm>
+#include <numeric>
#include <string>
#include <vector>
#include "gromacs/utility/stringutil.h"
#include "gromacs/utility/textwriter.h"
+namespace gmx
+{
+
+IndexGroupsAndNames::IndexGroupsAndNames(
+ const t_blocka &indexGroup, ArrayRef<char const * const> groupNames)
+ : indexGroup_ {indexGroup}
+{
+ std::copy(groupNames.begin(), groupNames.end(), std::back_inserter(groupNames_));
+ GMX_ASSERT(indexGroup_.nr == ssize(groupNames),
+ "Number of groups must match number of group names.");
+}
+
+bool IndexGroupsAndNames::containsGroupName(const std::string &groupName) const
+{
+ return std::any_of(std::begin(groupNames_), std::end(groupNames_),
+ [&groupName](const std::string &name){return equalCaseInsensitive(groupName, name); });
+}
+
+std::vector<index> IndexGroupsAndNames::indices(const std::string &groupName) const
+{
+ if (!containsGroupName(groupName))
+ {
+ GMX_THROW(InconsistentInputError(std::string("Group ") + groupName +
+ " referenced in the .mdp file was not found in the index file.\n"
+ "Group names must match either [moleculetype] names or custom index group\n"
+ "names, in which case you must supply an index file to the '-n' option\n"
+ "of grompp."));
+ }
+ const auto groupNamePosition = std::find_if(std::begin(groupNames_), std::end(groupNames_),
+ [&groupName](const std::string &name){return equalCaseInsensitive(groupName, name); });
+ const auto groupIndex = std::distance(std::begin(groupNames_), groupNamePosition);
+ const auto groupSize = indexGroup_.index[groupIndex+1] - indexGroup_.index[groupIndex];
+ std::vector<index> groupIndices(groupSize);
+ const auto startingIndex = indexGroup_.index[groupIndex];
+ std::iota(std::begin(groupIndices), std::end(groupIndices), startingIndex);
+ std::transform(std::begin(groupIndices), std::end(groupIndices), std::begin(groupIndices),
+ [blockLookup = indexGroup_.a](auto i){return blockLookup[i]; });
+ return groupIndices;
+}
+
+} // namespace gmx
+
/********************************************************************
* gmx_ana_indexgrps_t functions
********************************************************************/
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016,2018, by the GROMACS development team, led by
+ * Copyright (c) 2009,2010,2011,2012,2013,2014,2015,2016,2018,2019, 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 <cstdio>
#include <string>
+#include <vector>
#include "gromacs/topology/block.h"
+#include "gromacs/utility/arrayref.h"
namespace gmx
{
class TextWriter;
-}
+
+/*!\internal \brief Bundle index groups with their names.
+ *
+ * \note This class has to outlive the t_blocka &indexGroup it is given.
+ * We want to avoid a deep-copy due to the potentially very large data
+ * stored in the t_blocka.
+ *
+ * \todo update this class, once the two legacy data structures, t_blocka and
+ * the char ** of group names are refactored.
+ */
+class IndexGroupsAndNames
+{
+ public:
+
+ /*!\brief Construct from index group and group names
+ * \param[in] indexGroup
+ * \param[in] groupNames names of the index groups
+ */
+ IndexGroupsAndNames(const t_blocka &indexGroup, ArrayRef<char const *const> groupNames);
+
+ /*!\brief Return if a group name is contained in the groups.
+ *
+ * String comparison is case insensitive
+ *
+ * \param[in] groupName the group name to be queried
+ * \returns true if index group name is contained
+ */
+ bool containsGroupName(const std::string &groupName) const;
+
+ /*!\brief Return the integer indices of a group.
+ *
+ * If two index groups share a name, return the one found first.
+ *
+ * Indices may be empty.
+ *
+ * \param[in] groupName the name of the group whose indices shall be returned
+ * \returns atom indices of the selected index group
+ * \throws if groupName is not present as index group
+ */
+ std::vector<index> indices(const std::string &groupName) const;
+ private:
+ const t_blocka &indexGroup_;
+ std::vector<std::string> groupNames_;
+};
+
+} // namespace gmx
struct gmx_mtop_t;
#include "gromacs/selection/indexutil.h"
+#include <gmock/gmock.h>
#include <gtest/gtest.h>
#include "gromacs/topology/block.h"
#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/smalloc.h"
#include "testutils/refdata.h"
+#include "testutils/testasserts.h"
#include "toputils.h"
testUpdate(maxGroup, false, "Initialized");
}
+/***********************************************************************
+ * IndexGroupsAndNames tests
+ */
+
+class IndexGroupsAndNamesTest : public ::testing::Test
+{
+ public:
+ IndexGroupsAndNamesTest()
+ {
+ init_blocka(&blockA_);
+ addGroupToBlocka_(indicesGroupA_);
+ addGroupToBlocka_(indicesGroupB_);
+ addGroupToBlocka_(indicesGroupSecondA_);
+ addGroupToBlocka_(indicesGroupC_);
+
+ const char *const namesAsConstCharArray[4] = {
+ groupNames[0].c_str(), groupNames[1].c_str(),
+ groupNames[2].c_str(), groupNames[3].c_str()
+ };
+ indexGroupAndNames_ = std::make_unique<gmx::IndexGroupsAndNames>(blockA_, namesAsConstCharArray);
+ }
+ ~IndexGroupsAndNamesTest() override
+ {
+ done_blocka(&blockA_);
+ }
+ protected:
+ std::unique_ptr<gmx::IndexGroupsAndNames> indexGroupAndNames_;
+ const std::vector<std::string> groupNames = { "A", "B - Name", "A", "C" };
+ const std::vector<gmx::index> indicesGroupA_ = { };
+ const std::vector<gmx::index> indicesGroupB_ = { 1, 2 };
+ const std::vector<gmx::index> indicesGroupSecondA_ = { 5 };
+ const std::vector<gmx::index> indicesGroupC_ = { 10 };
+
+ private:
+ //! Add a new group to t_blocka
+ void addGroupToBlocka_(gmx::ArrayRef<const gmx::index> index)
+ {
+ srenew(blockA_.index, blockA_.nr + 2);
+ srenew(blockA_.a, blockA_.nra + index.size());
+ for (int i = 0; (i < index.ssize()); i++)
+ {
+ blockA_.a[blockA_.nra++] = index[i];
+ }
+ blockA_.nr++;
+ blockA_.index[blockA_.nr] = blockA_.nra;
+ }
+
+ t_blocka blockA_;
+};
+
+TEST_F(IndexGroupsAndNamesTest, containsNames)
+{
+ EXPECT_TRUE(indexGroupAndNames_->containsGroupName("a"));
+ EXPECT_TRUE(indexGroupAndNames_->containsGroupName("A"));
+ EXPECT_TRUE(indexGroupAndNames_->containsGroupName("B - Name"));
+ EXPECT_TRUE(indexGroupAndNames_->containsGroupName("b - Name"));
+ EXPECT_TRUE(indexGroupAndNames_->containsGroupName("B - naMe"));
+ EXPECT_TRUE(indexGroupAndNames_->containsGroupName("C"));
+ EXPECT_FALSE(indexGroupAndNames_->containsGroupName("D"));
+}
+
+TEST_F(IndexGroupsAndNamesTest, throwsWhenNameMissing)
+{
+ EXPECT_ANY_THROW(indexGroupAndNames_->indices("D"));
+}
+
+TEST_F(IndexGroupsAndNamesTest, groupIndicesCorrect)
+{
+ using ::testing::Pointwise;
+ using ::testing::Eq;
+ EXPECT_THAT(indicesGroupA_, Pointwise(Eq(), indexGroupAndNames_->indices("A")));
+ EXPECT_THAT(indicesGroupB_, Pointwise(Eq(), indexGroupAndNames_->indices("B - name")));
+ EXPECT_THAT(indicesGroupC_, Pointwise(Eq(), indexGroupAndNames_->indices("C")));
+}
+
+
} // namespace
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff