From: Christian Blau Date: Wed, 14 Aug 2019 10:30:40 +0000 (+0200) Subject: Add index group option to density fitting X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=8bec4630cd63c8db6b2899c9337f87f0a1ff943d;p=alexxy%2Fgromacs.git Add index group option to density fitting Adds the option to select a specific group for density fitting, stores the indices in the tpr file and restores them after reading. Adapted the densityfittingoptions test to test fixture to better reflect the increased required complexity for testing the storage of indices. refs #2282 Change-Id: I3e33508bbcadf968b0ddf978a73e7f6cf5753b5b --- diff --git a/src/gromacs/applied_forces/CMakeLists.txt b/src/gromacs/applied_forces/CMakeLists.txt index d5a61a6e07..d7ea182db6 100644 --- a/src/gromacs/applied_forces/CMakeLists.txt +++ b/src/gromacs/applied_forces/CMakeLists.txt @@ -38,6 +38,7 @@ gmx_add_libgromacs_sources( densityfittingforceprovider.cpp densityfittingoptions.cpp densityfittingoutputprovider.cpp + densityfittingparameters.cpp electricfield.cpp ) diff --git a/src/gromacs/applied_forces/densityfitting.cpp b/src/gromacs/applied_forces/densityfitting.cpp index 7b8b9fb2ea..de1d6adac0 100644 --- a/src/gromacs/applied_forces/densityfitting.cpp +++ b/src/gromacs/applied_forces/densityfitting.cpp @@ -45,11 +45,13 @@ #include "gromacs/mdrunutility/mdmodulenotification.h" #include "gromacs/mdtypes/imdmodule.h" +#include "gromacs/utility/keyvaluetreebuilder.h" #include "densityfittingforceprovider.h" #include "densityfittingoptions.h" #include "densityfittingoutputprovider.h" + namespace gmx { @@ -69,9 +71,42 @@ class DensityFitting final : public IMDModule { public: /*! \brief Construct the density fitting module. - * Allow the module to subscribe to notifications from MdModules + * + * \param[in] notifier allows the module to subscribe to notifications from MdModules. + * + * The density fitting code subscribes to these notifications: + * - setting atom group indices in the densityFittingOptions_ by + * taking a parmeter const IndexGroupsAndNames & + * - storing its internal parameters in a tpr file by writing to a + * key-value-tree during pre-processing by a function taking a + * KeyValueTreeObjectBuilder as parameter + * - reading its internal parameters from a key-value-tree during + * simulation setup by taking a const KeyValueTreeObject & parameter */ - explicit DensityFitting(MdModulesNotifier * /*notifier*/){} + explicit DensityFitting(MdModulesNotifier *notifier) + { + // Callbacks for several kinds of MdModuleNotification are created + // and subscribed, and will be dispatched correctly at run time + // based on the type of the parameter required by the lambda. + + // Setting atom group indices + const auto setFitGroupIndicesFunction = [this](const IndexGroupsAndNames &indexGroupsAndNames) { + densityFittingOptions_.setFitGroupIndices(indexGroupsAndNames); + }; + notifier->notifier_.subscribe(setFitGroupIndicesFunction); + + // Writing internal parameters during pre-processing + const auto writeInternalParametersFunction = [this](KeyValueTreeObjectBuilder treeBuilder) { + densityFittingOptions_.writeInternalParametersToKvt(treeBuilder); + }; + notifier->notifier_.subscribe(writeInternalParametersFunction); + + // Reading internal parameters during simulation setup + const auto readInternalParametersFunction = [this](const KeyValueTreeObject &tree) { + densityFittingOptions_.readInternalParametersFromKvt(tree); + }; + notifier->notifier_.subscribe(readInternalParametersFunction); + } //! From IMDModule; this class provides the mdpOptions itself IMdpOptionProvider *mdpOptionProvider() override { return &densityFittingOptions_; } diff --git a/src/gromacs/applied_forces/densityfittingoptions.cpp b/src/gromacs/applied_forces/densityfittingoptions.cpp index e0ce323e2d..7d5400d5b7 100644 --- a/src/gromacs/applied_forces/densityfittingoptions.cpp +++ b/src/gromacs/applied_forces/densityfittingoptions.cpp @@ -46,6 +46,8 @@ #include "gromacs/applied_forces/densityfitting.h" #include "gromacs/options/basicoptions.h" #include "gromacs/options/optionsection.h" +#include "gromacs/selection/indexutil.h" +#include "gromacs/utility/exceptions.h" #include "gromacs/utility/keyvaluetreebuilder.h" #include "gromacs/utility/keyvaluetreetransform.h" #include "gromacs/utility/strconvert.h" @@ -104,12 +106,20 @@ void addDensityFittingMdpOutputValue(KeyValueTreeObjectBuilder *builder, void DensityFittingOptions::initMdpTransform(IKeyValueTreeTransformRules * rules) { + const auto &stringIdentityTransform = [](std::string s){ + return s; + }; densityfittingMdpTransformFromString(rules, &fromStdString, c_activeTag_); + densityfittingMdpTransformFromString(rules, stringIdentityTransform, c_groupTag_); } void DensityFittingOptions::buildMdpOutput(KeyValueTreeObjectBuilder *builder) const { addDensityFittingMdpOutputValue(builder, parameters_.active_, c_activeTag_); + if (parameters_.active_) + { + addDensityFittingMdpOutputValue(builder, groupString_, c_groupTag_); + } } void DensityFittingOptions::initMdpOptions(IOptionsContainerWithSections *options) @@ -117,6 +127,7 @@ void DensityFittingOptions::initMdpOptions(IOptionsContainerWithSections *option auto section = options->addSection(OptionSection(DensityFittingModuleInfo::name_.c_str())); section.addOption(BooleanOption(c_activeTag_.c_str()).store(¶meters_.active_)); + section.addOption(StringOption(c_groupTag_.c_str()).store(&groupString_)); } bool DensityFittingOptions::active() const @@ -129,4 +140,40 @@ const DensityFittingParameters &DensityFittingOptions::buildParameters() return parameters_; } +void DensityFittingOptions::setFitGroupIndices(const IndexGroupsAndNames &indexGroupsAndNames) +{ + if (!parameters_.active_) + { + return; + } + parameters_.indices_ = indexGroupsAndNames.indices(groupString_); +} + +void DensityFittingOptions::writeInternalParametersToKvt(KeyValueTreeObjectBuilder treeBuilder) +{ + auto groupIndexAdder = treeBuilder.addUniformArray(DensityFittingModuleInfo::name_ + "-" + c_groupTag_); + for (const auto &indexValue : parameters_.indices_) + { + groupIndexAdder.addValue(indexValue); + } +} + +void DensityFittingOptions::readInternalParametersFromKvt(const KeyValueTreeObject &tree) +{ + if (!parameters_.active_) + { + return; + } + + if (!tree.keyExists(DensityFittingModuleInfo::name_ + "-" + c_groupTag_)) + { + GMX_THROW(InconsistentInputError( + "Cannot find atom index vector required for density guided simulation.")); + } + auto kvtIndexArray = tree[DensityFittingModuleInfo::name_ + "-" + c_groupTag_].asArray().values(); + parameters_.indices_.resize(kvtIndexArray.size()); + std::transform(std::begin(kvtIndexArray), std::end(kvtIndexArray), std::begin(parameters_.indices_), + [](const KeyValueTreeValue &val) { return val.cast(); }); +} + } // namespace gmx diff --git a/src/gromacs/applied_forces/densityfittingoptions.h b/src/gromacs/applied_forces/densityfittingoptions.h index 523c41f7c5..e0c7394269 100644 --- a/src/gromacs/applied_forces/densityfittingoptions.h +++ b/src/gromacs/applied_forces/densityfittingoptions.h @@ -51,6 +51,10 @@ namespace gmx { +class IndexGroupsAndNames; +class KeyValueTreeObject; +class KeyValueTreeBuilder; + /*! \internal * \brief Input data storage for density fitting */ @@ -80,8 +84,29 @@ class DensityFittingOptions final : public IMdpOptionProvider //! Process input options to parameters, including input file reading. const DensityFittingParameters &buildParameters(); + /*! \brief Evaluate and store atom indices. + * + * During pre-processing, use the group string from the options to + * evaluate the indices of the atoms to be subject to forces from this + * module. + */ + void setFitGroupIndices(const IndexGroupsAndNames &indexGroupsAndNames); + + //! Store the paramers that are not mdp options in the tpr file + void writeInternalParametersToKvt(KeyValueTreeObjectBuilder treeBuilder); + + //! Set the internal parameters that are stored in the tpr file + void readInternalParametersFromKvt(const KeyValueTreeObject &tree); + private: const std::string c_activeTag_ = "active"; + + /*! \brief Denote the .mdp option that defines the group of fit atoms. + * \note Changing this string will break .tpr backwards compability + */ + const std::string c_groupTag_ = "group"; + std::string groupString_ = "protein"; + DensityFittingParameters parameters_; }; diff --git a/src/gromacs/applied_forces/densityfittingparameters.cpp b/src/gromacs/applied_forces/densityfittingparameters.cpp new file mode 100644 index 0000000000..f390d439e1 --- /dev/null +++ b/src/gromacs/applied_forces/densityfittingparameters.cpp @@ -0,0 +1,60 @@ +/* + * This file is part of the GROMACS molecular simulation package. + * + * Copyright (c) 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. + * + * 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. + */ + +#include "gmxpre.h" + +#include "densityfittingparameters.h" + +namespace gmx +{ +bool operator==(const DensityFittingParameters &lhs, const DensityFittingParameters &rhs) +{ + if (lhs.active_ != rhs.active_) + { + return false; + } + if (lhs.indices_ != rhs.indices_) + { + return false; + } + return true; +} + +bool operator!=(const DensityFittingParameters &lhs, const DensityFittingParameters &rhs) +{ + return !(lhs == rhs); +} + +} // namespace gmx diff --git a/src/gromacs/applied_forces/densityfittingparameters.h b/src/gromacs/applied_forces/densityfittingparameters.h index e2df371291..25fcc54b56 100644 --- a/src/gromacs/applied_forces/densityfittingparameters.h +++ b/src/gromacs/applied_forces/densityfittingparameters.h @@ -42,6 +42,10 @@ #ifndef GMX_APPLIED_FORCES_DENSITYFITTINGPARAMETERS_H #define GMX_APPLIED_FORCES_DENSITYFITTINGPARAMETERS_H +#include + +#include "gromacs/utility/basedefinitions.h" + namespace gmx { @@ -53,9 +57,27 @@ namespace gmx struct DensityFittingParameters { //! Indicate if density fitting is active - bool active_ = false; + bool active_ = false; + //! Indices of the atoms that shall be fit to the density + std::vector indices_; }; +/*! \brief Check if two structs holding density fitting parameters are equal. + * + * \param[in] lhs left hand side to be compared + * \param[in] rhs right hand side to be compared + * \returns true if all elements in DensityFittingParameters are equal, else false + */ +bool operator==(const DensityFittingParameters &lhs, const DensityFittingParameters &rhs); + +/*! \brief Check if two structs holding density fitting parameters are not equal. + * + * \param[in] lhs left hand side to be compared + * \param[in] rhs right hand side to be compared + * \returns true if lhs is not equal rhs + */ +bool operator!=(const DensityFittingParameters &lhs, const DensityFittingParameters &rhs); + } // namespace gmx #endif // GMX_APPLIED_FORCES_DENSITYFITTINGPARAMETERS_H diff --git a/src/gromacs/applied_forces/tests/densityfittingoptions.cpp b/src/gromacs/applied_forces/tests/densityfittingoptions.cpp index 358a9f42a2..a6c7f57d1d 100644 --- a/src/gromacs/applied_forces/tests/densityfittingoptions.cpp +++ b/src/gromacs/applied_forces/tests/densityfittingoptions.cpp @@ -43,13 +43,18 @@ #include "gromacs/applied_forces/densityfittingoptions.h" +#include +#include + #include #include "gromacs/options/options.h" #include "gromacs/options/treesupport.h" +#include "gromacs/selection/indexutil.h" #include "gromacs/utility/keyvaluetreebuilder.h" #include "gromacs/utility/keyvaluetreemdpwriter.h" #include "gromacs/utility/keyvaluetreetransform.h" +#include "gromacs/utility/smalloc.h" #include "gromacs/utility/stringcompare.h" #include "gromacs/utility/stringstream.h" #include "gromacs/utility/textwriter.h" @@ -57,48 +62,92 @@ #include "testutils/testasserts.h" #include "testutils/testmatchers.h" - namespace gmx { namespace { -TEST(DensityFittingOptionsTest, DefaultParameters) +class DensityFittingOptionsTest : public ::testing::Test { - DensityFittingOptions densityFittingOptions; - EXPECT_FALSE(densityFittingOptions.buildParameters().active_); + public: + DensityFittingOptionsTest() + { + init_blocka(&defaultGroups_); + } + ~DensityFittingOptionsTest() override + { + done_blocka(&defaultGroups_); + } + + void setFromMdpValues(const KeyValueTreeObject &densityFittingMdpValues) + { + // set up options + Options densityFittingModuleOptions; + densityFittingOptions_.initMdpOptions(&densityFittingModuleOptions); + + // Add rules to transform mdp inputs to densityFittingModule data + KeyValueTreeTransformer transform; + transform.rules()->addRule().keyMatchType("/", StringCompareType::CaseAndDashInsensitive); + + densityFittingOptions_.initMdpTransform(transform.rules()); + + // Execute the transform on the mdpValues + auto transformedMdpValues = transform.transform(densityFittingMdpValues, nullptr); + assignOptionsFromKeyValueTree(&densityFittingModuleOptions, transformedMdpValues.object(), nullptr); + } + + KeyValueTreeObject densityFittingSetActiveAsMdpValues() + { + // Prepare MDP inputs + KeyValueTreeBuilder mdpValueBuilder; + mdpValueBuilder.rootObject().addValue("density-guided-simulation-active", + std::string("yes")); + return mdpValueBuilder.build(); + } + + IndexGroupsAndNames genericIndexGroupsAndNames() + { + done_blocka(&defaultGroups_); + stupid_fill_blocka(&defaultGroups_, 3); + std::vector groupNames = { "A", "protein", "C" }; + const char *const namesAsConstChar[3] + = { groupNames[0].c_str(), groupNames[1].c_str(), groupNames[2].c_str() }; + return {defaultGroups_, namesAsConstChar}; + } + + IndexGroupsAndNames differingIndexGroupsAndNames() + { + done_blocka(&defaultGroups_); + stupid_fill_blocka(&defaultGroups_, 3); + std::vector groupNames = { "protein", "C", "A"}; + const char *const namesAsConstChar[3] + = { groupNames[0].c_str(), groupNames[1].c_str(), groupNames[2].c_str() }; + return { defaultGroups_, namesAsConstChar }; + } + + void mangleInternalParameters() + { + densityFittingOptions_.setFitGroupIndices(differingIndexGroupsAndNames()); + } + protected: + t_blocka defaultGroups_; + DensityFittingOptions densityFittingOptions_; +}; + +TEST_F(DensityFittingOptionsTest, DefaultParameters) +{ + EXPECT_FALSE(densityFittingOptions_.buildParameters().active_); } -TEST(DensityFittingOptionsTest, OptionSetsActive) +TEST_F(DensityFittingOptionsTest, OptionSetsActive) { - DensityFittingOptions densityFittingOptions; - EXPECT_FALSE(densityFittingOptions.buildParameters().active_); - - // Prepare MDP inputs - KeyValueTreeBuilder mdpValueBuilder; - mdpValueBuilder.rootObject().addValue("density-guided-simulation-active", - std::string("yes")); - KeyValueTreeObject densityFittingMdpValues = mdpValueBuilder.build(); - - // set up options - Options densityFittingModuleOptions; - densityFittingOptions.initMdpOptions(&densityFittingModuleOptions); - - // Add rules to transform mdp inputs to densityFittingModule data - KeyValueTreeTransformer transform; - transform.rules()->addRule().keyMatchType("/", StringCompareType::CaseAndDashInsensitive); - - densityFittingOptions.initMdpTransform(transform.rules()); - - // Execute the transform on the mdpValues - auto transformedMdpValues = transform.transform(densityFittingMdpValues, nullptr); - assignOptionsFromKeyValueTree(&densityFittingModuleOptions, transformedMdpValues.object(), nullptr); - - EXPECT_TRUE(densityFittingOptions.buildParameters().active_); + EXPECT_FALSE(densityFittingOptions_.buildParameters().active_); + setFromMdpValues(densityFittingSetActiveAsMdpValues()); + EXPECT_TRUE(densityFittingOptions_.buildParameters().active_); } -TEST(DensityFittingOptionsTest, OutputDefaultValues) +TEST_F(DensityFittingOptionsTest, OutputDefaultValues) { // Transform module data into a flat key-value tree for output. @@ -106,9 +155,7 @@ TEST(DensityFittingOptionsTest, OutputDefaultValues) KeyValueTreeBuilder builder; KeyValueTreeObjectBuilder builderObject = builder.rootObject(); - DensityFittingOptions densityFittingOptions; - densityFittingOptions.buildMdpOutput(&builderObject); - + densityFittingOptions_.buildMdpOutput(&builderObject); { TextWriter writer(&stream); writeKeyValueTreeAsMdp(&writer, builder.build()); @@ -118,6 +165,73 @@ TEST(DensityFittingOptionsTest, OutputDefaultValues) EXPECT_EQ(stream.toString(), std::string("density-guided-simulation-active = false\n")); } +TEST_F(DensityFittingOptionsTest, CanConvertGroupStringToIndexGroup) +{ + setFromMdpValues(densityFittingSetActiveAsMdpValues()); + + const auto indexGroupAndNames = genericIndexGroupsAndNames(); + densityFittingOptions_.setFitGroupIndices(indexGroupAndNames); + + EXPECT_EQ(1, densityFittingOptions_.buildParameters().indices_.size()); + EXPECT_EQ(1, densityFittingOptions_.buildParameters().indices_[0]); +} + +TEST_F(DensityFittingOptionsTest, InternalsToKvt) +{ + // stores the default internal options + DensityFittingOptions densityFittingOptions; + KeyValueTreeBuilder builder; + densityFittingOptions.writeInternalParametersToKvt(builder.rootObject()); + const auto kvtTree = builder.build(); + EXPECT_TRUE(kvtTree.keyExists("density-guided-simulation-group")); + EXPECT_TRUE(kvtTree["density-guided-simulation-group"].isArray()); + auto storedIndex = kvtTree["density-guided-simulation-group"].asArray().values(); + + EXPECT_EQ(0, storedIndex.size()); +} + +TEST_F(DensityFittingOptionsTest, KvtToInternal) +{ + setFromMdpValues(densityFittingSetActiveAsMdpValues()); + + KeyValueTreeBuilder builder; + auto addedArray = builder.rootObject().addUniformArray("density-guided-simulation-group"); + addedArray.addValue(1); + addedArray.addValue(15); + const auto tree = builder.build(); + + densityFittingOptions_.readInternalParametersFromKvt(tree); + + EXPECT_EQ(2, densityFittingOptions_.buildParameters().indices_.size()); + EXPECT_EQ(1, densityFittingOptions_.buildParameters().indices_[0]); + EXPECT_EQ(15, densityFittingOptions_.buildParameters().indices_[1]); +} + +TEST_F(DensityFittingOptionsTest, RoundTripForInternalsIsIdempotent) +{ + setFromMdpValues(densityFittingSetActiveAsMdpValues()); + { + const IndexGroupsAndNames indexGroupAndNames = genericIndexGroupsAndNames(); + densityFittingOptions_.setFitGroupIndices(indexGroupAndNames); + } + + DensityFittingParameters parametersBefore = densityFittingOptions_.buildParameters(); + + KeyValueTreeBuilder builder; + densityFittingOptions_.writeInternalParametersToKvt(builder.rootObject()); + const auto inputTree = builder.build(); + + mangleInternalParameters(); + + DensityFittingParameters parametersAfter = densityFittingOptions_.buildParameters(); + EXPECT_NE(parametersBefore, parametersAfter); + + densityFittingOptions_.readInternalParametersFromKvt(inputTree); + + parametersAfter = densityFittingOptions_.buildParameters(); + EXPECT_EQ(parametersBefore, parametersAfter); +} + } // namespace } // namespace gmx diff --git a/src/gromacs/gmxpreprocess/grompp.cpp b/src/gromacs/gmxpreprocess/grompp.cpp index f15c17ba7d..6cedcc3441 100644 --- a/src/gromacs/gmxpreprocess/grompp.cpp +++ b/src/gromacs/gmxpreprocess/grompp.cpp @@ -2425,7 +2425,7 @@ int gmx_grompp(int argc, char *argv[]) { gmx::KeyValueTreeBuilder internalParameterBuilder; - mdModules.notifier().notifier_.notify(&internalParameterBuilder); + mdModules.notifier().notifier_.notify(internalParameterBuilder.rootObject()); ir->internalParameters = std::make_unique(internalParameterBuilder.build()); } diff --git a/src/gromacs/mdrun/mdmodules.cpp b/src/gromacs/mdrun/mdmodules.cpp index 2f1119c06e..1e1104b2e6 100644 --- a/src/gromacs/mdrun/mdmodules.cpp +++ b/src/gromacs/mdrun/mdmodules.cpp @@ -93,6 +93,13 @@ class MDModules::Impl : public IMDOutputProvider densityFitting_->outputProvider()->finishOutput(); } + /*! \brief Manages callbacks and notifies the MD modules. + * + * \note The notifier must be constructed before the modules and shall + * not be destructed before the modules are destructed. + */ + MdModulesNotifier notifier_; + std::unique_ptr densityFitting_; std::unique_ptr field_; std::unique_ptr forceProviders_; @@ -109,9 +116,6 @@ class MDModules::Impl : public IMDOutputProvider * \todo include field_ in modules_ */ std::vector< std::shared_ptr > modules_; - - //! Manages resources and notifies the MD modules when available - MdModulesNotifier notifier_; }; MDModules::MDModules() : impl_(new Impl) diff --git a/src/gromacs/mdrunutility/mdmodulenotification.h b/src/gromacs/mdrunutility/mdmodulenotification.h index dea6e1637f..d7851fdd5b 100644 --- a/src/gromacs/mdrunutility/mdmodulenotification.h +++ b/src/gromacs/mdrunutility/mdmodulenotification.h @@ -195,7 +195,7 @@ struct registerMdModuleNotification }; class KeyValueTreeObject; -class KeyValueTreeBuilder; +class KeyValueTreeObjectBuilder; class LocalAtomSetManager; class IndexGroupsAndNames; @@ -205,7 +205,7 @@ struct MdModulesNotifier registerMdModuleNotification< const t_commrec &, IndexGroupsAndNames, - KeyValueTreeBuilder *, + KeyValueTreeObjectBuilder, const KeyValueTreeObject &, LocalAtomSetManager *>::type notifier_; };