Add reference density file option to densityfitting

Add option to density fitting to read reference density for
density fitting code.

refs #2282

Change-Id: I598b679a3c1c63b913a104d7ccff8e1643b3946d

diff --git a/src/gromacs/applied_forces/densityfitting.cpp b/src/gromacs/applied_forces/densityfitting.cpp
index 357e8ddfa610fb739d20a455b54d4451ab4ebabb..a276e35fcd8d170c88a6d176bf134d276441e85e 100644 (file)
--- a/src/gromacs/applied_forces/densityfitting.cpp
+++ b/src/gromacs/applied_forces/densityfitting.cpp
@@ -46,6 +46,8 @@
 #include <memory>
 
 #include "gromacs/domdec/localatomsetmanager.h"
+#include "gromacs/fileio/mrcdensitymap.h"
+#include "gromacs/math/multidimarray.h"
 #include "gromacs/mdtypes/imdmodule.h"
 #include "gromacs/utility/classhelpers.h"
 #include "gromacs/utility/exceptions.h"
@@ -56,7 +58,6 @@
 #include "densityfittingoptions.h"
 #include "densityfittingoutputprovider.h"
 
-
 namespace gmx
 {
 
@@ -84,6 +85,7 @@ class DensityFittingSimulationParameterSetup
 {
     public:
         DensityFittingSimulationParameterSetup() = default;
+
         /*! \brief Set the local atom set for the density fitting.
          * \param[in] localAtomSet of atoms to be fitted
          */
@@ -101,7 +103,7 @@ class DensityFittingSimulationParameterSetup
             if (localAtomSet_ == nullptr)
             {
                 GMX_THROW(
-                        InternalError("Transformation to reference density not set for density "
+                        InternalError("Local atom set is not set for density "
                                       "guided simulation."));
             }
             return *localAtomSet_;
@@ -124,22 +126,40 @@ class DensityFittingSimulationParameterSetup
          * \throws InternalError if reference density is not set
          * \returns the reference density
          */
-        const basic_mdspan<const float, dynamicExtents3D> &referenceDensity() const
+        basic_mdspan<const float, dynamicExtents3D>
+        referenceDensity() const
         {
             if (referenceDensity_ == nullptr)
             {
                 GMX_THROW(InternalError("Reference density not set for density guided simulation."));
             }
-            return *referenceDensity_;
+            return referenceDensity_->asConstView();
+        }
+
+        /*! \brief Reads the reference density from file.
+         *
+         * Reads and check file, then set and communicate the internal
+         * parameters related to the reference density with the file data.
+         *
+         * \throws FileIOError if reading from file was not successful
+         */
+        void readReferenceDensityFromFile(const std::string &referenceDensityFileName)
+        {
+            MrcDensityMapOfFloatFromFileReader reader(referenceDensityFileName);
+            referenceDensity_ = std::make_unique<MultiDimArray<std::vector<float>, dynamicExtents3D> >
+                    (reader.densityDataCopy());
+            transformationToDensityLattice_
+                = std::make_unique<TranslateAndScale>(reader.transformationToDensityLattice());
         }
 
     private:
         //! The reference density to fit to
-        std::unique_ptr < basic_mdspan < const float, dynamicExtents3D>> referenceDensity_;
+        std::unique_ptr<MultiDimArray<std::vector<float>, dynamicExtents3D> > referenceDensity_;
         //! The coordinate transformation into the reference density
         std::unique_ptr<TranslateAndScale> transformationToDensityLattice_;
         //! The local atom set to act on
         std::unique_ptr<LocalAtomSet>      localAtomSet_;
+
         GMX_DISALLOW_COPY_AND_ASSIGN(DensityFittingSimulationParameterSetup);
 };
 
@@ -198,7 +218,7 @@ class DensityFitting final : public IMDModule
 
         }
 
-        //! From IMDModule; this class provides the mdpOptions itself
+        //! From IMDModule
         IMdpOptionProvider *mdpOptionProvider() override { return &densityFittingOptions_; }
 
         //! Add this module to the force providers if active
@@ -207,6 +227,7 @@ class DensityFitting final : public IMDModule
             if (densityFittingOptions_.active())
             {
                 const auto &parameters = densityFittingOptions_.buildParameters();
+                densityFittingSimulationParameters_.readReferenceDensityFromFile(densityFittingOptions_.referenceDensityFileName());
                 forceProvider_ = std::make_unique<DensityFittingForceProvider>(
                             parameters,
                             densityFittingSimulationParameters_.referenceDensity(),
@@ -243,6 +264,7 @@ class DensityFitting final : public IMDModule
          * simulation setup time.
          */
         DensityFittingSimulationParameterSetup       densityFittingSimulationParameters_;
+
         GMX_DISALLOW_COPY_AND_ASSIGN(DensityFitting);
 };
 

diff --git a/src/gromacs/applied_forces/densityfittingoptions.cpp b/src/gromacs/applied_forces/densityfittingoptions.cpp
index 68f1fb742b7eb4aed0672ef6478021ccfe2c1cbc..ccdb0043de503bc6a878d130f0fe527efc34efa8 100644 (file)
--- a/src/gromacs/applied_forces/densityfittingoptions.cpp
+++ b/src/gromacs/applied_forces/densityfittingoptions.cpp
@@ -119,6 +119,7 @@ void DensityFittingOptions::initMdpTransform(IKeyValueTreeTransformRules * rules
     densityfittingMdpTransformFromString<real>(rules, &fromStdString<real>, c_forceConstantTag_);
     densityfittingMdpTransformFromString<real>(rules, &fromStdString<real>, c_gaussianTransformSpreadingWidthTag_);
     densityfittingMdpTransformFromString<real>(rules, &fromStdString<real>, c_gaussianTransformSpreadingRangeInMultiplesOfWidthTag_);
+    densityfittingMdpTransformFromString<std::string>(rules, stringIdentityTransform, c_referenceDensityFileNameTag_);
 }
 
 void DensityFittingOptions::buildMdpOutput(KeyValueTreeObjectBuilder *builder) const
@@ -138,6 +139,7 @@ void DensityFittingOptions::buildMdpOutput(KeyValueTreeObjectBuilder *builder) c
         addDensityFittingMdpOutputValue(builder, parameters_.forceConstant_, c_forceConstantTag_);
         addDensityFittingMdpOutputValue(builder, parameters_.gaussianTransformSpreadingWidth_, c_gaussianTransformSpreadingWidthTag_);
         addDensityFittingMdpOutputValue(builder, parameters_.gaussianTransformSpreadingRangeInMultiplesOfWidth_, c_gaussianTransformSpreadingRangeInMultiplesOfWidthTag_);
+        addDensityFittingMdpOutputValue(builder, referenceDensityFileName_, c_referenceDensityFileNameTag_);
     }
 }
 
@@ -159,6 +161,7 @@ void DensityFittingOptions::initMdpOptions(IOptionsContainerWithSections *option
     section.addOption(RealOption(c_forceConstantTag_.c_str()).store(&parameters_.forceConstant_));
     section.addOption(RealOption(c_gaussianTransformSpreadingWidthTag_.c_str()).store(&parameters_.gaussianTransformSpreadingWidth_));
     section.addOption(RealOption(c_gaussianTransformSpreadingRangeInMultiplesOfWidthTag_.c_str()).store(&parameters_.gaussianTransformSpreadingRangeInMultiplesOfWidth_));
+    section.addOption(StringOption(c_referenceDensityFileNameTag_.c_str()).store(&referenceDensityFileName_));
 }
 
 bool DensityFittingOptions::active() const
@@ -207,4 +210,8 @@ void DensityFittingOptions::readInternalParametersFromKvt(const KeyValueTreeObje
                    [](const KeyValueTreeValue &val) { return val.cast<index>(); });
 }
 
+const std::string &DensityFittingOptions::referenceDensityFileName() const
+{
+    return referenceDensityFileName_;
+}
 } // namespace gmx

diff --git a/src/gromacs/applied_forces/densityfittingoptions.h b/src/gromacs/applied_forces/densityfittingoptions.h
index d57ab61b5e76156dfea26c90ae8724e38e96d4d9..fa5ed45e7c89bea4c037bacb05c98212b18cb782 100644 (file)
--- a/src/gromacs/applied_forces/densityfittingoptions.h
+++ b/src/gromacs/applied_forces/densityfittingoptions.h
@@ -98,25 +98,32 @@ class DensityFittingOptions final : public IMdpOptionProvider
         //! Set the internal parameters that are stored in the tpr file
         void readInternalParametersFromKvt(const KeyValueTreeObject &tree);
 
+        //! Return the file name of the reference density
+        const std::string &referenceDensityFileName() const;
+
     private:
-        const std::string        c_activeTag_ = "active";
+        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";
+        const std::string c_groupTag_  = "group";
+        std::string       groupString_ = "protein";
 
-        const std::string        c_similarityMeasureTag_ = "similarity-measure";
+        const std::string c_similarityMeasureTag_ = "similarity-measure";
 
-        const std::string        c_amplitudeMethodTag_ = "amplitude-method";
+        const std::string c_amplitudeMethodTag_ = "amplitude-method";
 
-        const std::string        c_forceConstantTag_ = "force-constant";
+        const std::string c_forceConstantTag_ = "force-constant";
 
-        const std::string        c_gaussianTransformSpreadingWidthTag_ = "gaussian-transform-spreading-width";
-        const std::string        c_gaussianTransformSpreadingRangeInMultiplesOfWidthTag_
+        const std::string c_gaussianTransformSpreadingWidthTag_ = "gaussian-transform-spreading-width";
+        const std::string c_gaussianTransformSpreadingRangeInMultiplesOfWidthTag_
             = "gaussian-transform-spreading-range-in-multiples-of-width";
 
+        const std::string c_referenceDensityFileNameTag_ = "reference-density-filename";
+        std::string       referenceDensityFileName_      = "reference.mrc";
+
+
         DensityFittingParameters parameters_;
 };
 

diff --git a/src/gromacs/applied_forces/tests/densityfitting.cpp b/src/gromacs/applied_forces/tests/densityfitting.cpp
index 0bbe272400ce8dffac688bee6a66ce885ba9e433..c0a188fd0bb49ba2743edc4f1c61e95c289984f7 100644 (file)
--- a/src/gromacs/applied_forces/tests/densityfitting.cpp
+++ b/src/gromacs/applied_forces/tests/densityfitting.cpp
@@ -64,6 +64,7 @@
 #include "gromacs/utility/stringcompare.h"
 
 #include "testutils/testasserts.h"
+#include "testutils/testfilemanager.h"
 #include "testutils/testmatchers.h"
 
 namespace gmx
@@ -72,35 +73,79 @@ namespace gmx
 namespace
 {
 
-TEST(DensityFittingTest, ForceOnSingleOption)
+class DensityFittingTest : public ::testing::Test
 {
-    MdModulesNotifier notifier;
-    auto densityFittingModule(DensityFittingModuleInfo::create(&notifier));
+    public:
 
-    // Prepare MDP inputs
-    KeyValueTreeBuilder      mdpValueBuilder;
-    mdpValueBuilder.rootObject().addValue("density-guided-simulation-active", std::string("yes"));
-    KeyValueTreeObject       densityFittingMdpValues = mdpValueBuilder.build();
+        void addMdpOptionDensityFittingActive()
+        {
+            mdpValueBuilder_.rootObject().addValue("density-guided-simulation-active", std::string("yes"));
+        }
+
+        void addMdpOptionReferenceDensity()
+        {
+            mdpValueBuilder_.rootObject().addValue(
+                    "density-guided-simulation-reference-density-filename",
+                    std::string(test::TestFileManager::getInputFilePath("ellipsoid-density.mrc")));
+        }
+
+        //! build an mdp options tree that sets the options for the density fitting module
+        void makeDensityFittingModuleWithSetOptions()
+        {
+            KeyValueTreeObject mdpOptionsTree = mdpValueBuilder_.build();
+
+            densityFittingModule_ = DensityFittingModuleInfo::create(&notifier_);
+
+            // set up options
+            Options densityFittingModuleOptions;
+            densityFittingModule_->mdpOptionProvider()->initMdpOptions(&densityFittingModuleOptions);
 
-    // set up options
-    Options densityFittingModuleOptions;
-    densityFittingModule->mdpOptionProvider()->initMdpOptions(&densityFittingModuleOptions);
+            // Add rules to transform mdp inputs to densityFittingModule data
+            KeyValueTreeTransformer transform;
+            transform.rules()->addRule().keyMatchType("/", StringCompareType::CaseAndDashInsensitive);
+            densityFittingModule_->mdpOptionProvider()->initMdpTransform(transform.rules());
 
-    // Add rules to transform mdp inputs to densityFittingModule data
-    KeyValueTreeTransformer transform;
-    transform.rules()->addRule().keyMatchType("/", StringCompareType::CaseAndDashInsensitive);
-    densityFittingModule->mdpOptionProvider()->initMdpTransform(transform.rules());
+            // Execute the transform on the mdpValues
+            auto transformedMdpValues = transform.transform(mdpOptionsTree, nullptr);
+            assignOptionsFromKeyValueTree(&densityFittingModuleOptions, transformedMdpValues.object(), nullptr);
+        }
 
-    // Execute the transform on the mdpValues
-    auto transformedMdpValues = transform.transform(densityFittingMdpValues, nullptr);
-    assignOptionsFromKeyValueTree(&densityFittingModuleOptions, transformedMdpValues.object(), nullptr);
+        void intializeForceProviders()
+        {
+            densityFittingModule_->initForceProviders(&densityFittingForces_);
+        }
+
+    private:
+
+        KeyValueTreeBuilder        mdpValueBuilder_;
+        MdModulesNotifier          notifier_;
+        ForceProviders             densityFittingForces_;
+        std::unique_ptr<IMDModule> densityFittingModule_;
+};
+
+TEST_F(DensityFittingTest, ForceProviderLackingInputThrows)
+{
+    // Prepare MDP inputs
+    addMdpOptionDensityFittingActive();
+
+    makeDensityFittingModuleWithSetOptions();
 
     // Build the force provider, once all input data is gathered
-    ForceProviders densityFittingForces;
-    EXPECT_ANY_THROW(densityFittingModule->initForceProviders(&densityFittingForces));
+    EXPECT_ANY_THROW(intializeForceProviders());
+}
+
+TEST_F(DensityFittingTest, SingleAtom)
+{
+    // Prepare MDP inputs
+    addMdpOptionDensityFittingActive();
+    addMdpOptionReferenceDensity();
 
+    makeDensityFittingModuleWithSetOptions();
+
+    EXPECT_ANY_THROW(intializeForceProviders());
 }
 
+
 } // namespace
 
 } // namespace gmx

diff --git a/src/gromacs/applied_forces/tests/densityfittingoptions.cpp b/src/gromacs/applied_forces/tests/densityfittingoptions.cpp
index 7acd948da0ff2fd657046205f1cd7b3902ec6292..bd5d6d722d962c6880e78fd05a331d22c2c03f37 100644 (file)
--- a/src/gromacs/applied_forces/tests/densityfittingoptions.cpp
+++ b/src/gromacs/applied_forces/tests/densityfittingoptions.cpp
@@ -196,6 +196,7 @@ TEST_F(DensityFittingOptionsTest, OutputDefaultValuesWhenActive)
         "density-guided-simulation-force-constant = 1e+09\n"
         "density-guided-simulation-gaussian-transform-spreading-width = 0.2\n"
         "density-guided-simulation-gaussian-transform-spreading-range-in-multiples-of-width = 4\n"
+        "density-guided-simulation-reference-density-filename = reference.mrc\n"
         };
 
     EXPECT_EQ(expected, stream.toString());

diff --git a/src/gromacs/mdrun/runner.cpp b/src/gromacs/mdrun/runner.cpp
index 717c8c578c248167b5ad9c74e93e23ebe3b22c16..3ba528da84c3a7d2559f7a641caf69577c0c2406 100644 (file)
--- a/src/gromacs/mdrun/runner.cpp
+++ b/src/gromacs/mdrun/runner.cpp
@@ -992,7 +992,6 @@ int Mdrunner::mdrunner()
                                            &mtop, inputrec,
                                            box, positionsFromStatePointer(globalState.get()),
                                            &atomSets);
-        mdModulesNotifier.notify(&atomSets);
         // Note that local state still does not exist yet.
     }
     else
@@ -1279,6 +1278,7 @@ int Mdrunner::mdrunner()
     if (thisRankHasDuty(cr, DUTY_PP))
     {
         mdModulesNotifier.notify(*cr);
+        mdModulesNotifier.notify(&atomSets);
         /* Initiate forcerecord */
         fr                 = new t_forcerec;
         fr->forceProviders = mdModules_->initForceProviders();