src/gromacs/applied_forces/densityfitting.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Declares data structure and utilities for density fitting
  38  *
  39  * \author Christian Blau <blau@kth.se>
  40  * \ingroup module_applied_forces
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "densityfitting.h"
  45
  46 #include <memory>
  47
  48 #include "gromacs/domdec/localatomsetmanager.h"
  49 #include "gromacs/fileio/mrcdensitymap.h"
  50 #include "gromacs/math/multidimarray.h"
  51 #include "gromacs/mdtypes/imdmodule.h"
  52 #include "gromacs/utility/classhelpers.h"
  53 #include "gromacs/utility/exceptions.h"
  54 #include "gromacs/utility/keyvaluetreebuilder.h"
  55 #include "gromacs/utility/mdmodulenotification.h"
  56
  57 #include "densityfittingforceprovider.h"
  58 #include "densityfittingoptions.h"
  59 #include "densityfittingoutputprovider.h"
  60
  61 namespace gmx
  62 {
  63
  64 class IMdpOptionProvider;
  65 class DensityFittingForceProvider;
  66
  67 namespace
  68 {
  69
  70 /*! \internal
  71  * \brief Collect density fitting parameters only available during simulation setup.
  72  *
  73  * \todo Implement builder pattern that will not use unqiue_ptr to check if
  74  *       parameters have been set or not.
  75  *
  76  * To build the density fitting force provider during simulation setup,
  77  * the DensityFitting class needs access to parameters that become available
  78  * only during simulation setup.
  79  *
  80  * This class collects these parameters via MdModuleNotifications in the
  81  * simulation setup phase and provides a check if all necessary parameters have
  82  * been provided.
  83  */
  84 class DensityFittingSimulationParameterSetup
  85 {
  86     public:
  87         DensityFittingSimulationParameterSetup() = default;
  88
  89         /*! \brief Set the local atom set for the density fitting.
  90          * \param[in] localAtomSet of atoms to be fitted
  91          */
  92         void setLocalAtomSet(const LocalAtomSet &localAtomSet)
  93         {
  94             localAtomSet_ = std::make_unique<LocalAtomSet>(localAtomSet);
  95         }
  96
  97         /*! \brief Return local atom set for density fitting.
  98          * \throws InternalError if local atom set is not set
  99          * \returns local atom set for density fitting
 100          */
 101         const LocalAtomSet &localAtomSet() const
 102         {
 103             if (localAtomSet_ == nullptr)
 104             {
 105                 GMX_THROW(
 106                         InternalError("Local atom set is not set for density "
 107                                       "guided simulation."));
 108             }
 109             return *localAtomSet_;
 110         }
 111
 112         /*! \brief Return transformation into density lattice.
 113          * \throws InternalError if transformation into density lattice is not set
 114          * \returns transformation into density lattice
 115          */
 116         const TranslateAndScale &transformationToDensityLattice() const
 117         {
 118             if (transformationToDensityLattice_ == nullptr)
 119             {
 120                 GMX_THROW(
 121                         InternalError("Transformation to reference density not set for density guided simulation."));
 122             }
 123             return *transformationToDensityLattice_;
 124         }
 125         /*! \brief Return reference density
 126          * \throws InternalError if reference density is not set
 127          * \returns the reference density
 128          */
 129         basic_mdspan<const float, dynamicExtents3D>
 130         referenceDensity() const
 131         {
 132             if (referenceDensity_ == nullptr)
 133             {
 134                 GMX_THROW(InternalError("Reference density not set for density guided simulation."));
 135             }
 136             return referenceDensity_->asConstView();
 137         }
 138
 139         /*! \brief Reads the reference density from file.
 140          *
 141          * Reads and check file, then set and communicate the internal
 142          * parameters related to the reference density with the file data.
 143          *
 144          * \throws FileIOError if reading from file was not successful
 145          */
 146         void readReferenceDensityFromFile(const std::string &referenceDensityFileName)
 147         {
 148             MrcDensityMapOfFloatFromFileReader reader(referenceDensityFileName);
 149             referenceDensity_ = std::make_unique<MultiDimArray<std::vector<float>, dynamicExtents3D> >
 150                     (reader.densityDataCopy());
 151             transformationToDensityLattice_
 152                 = std::make_unique<TranslateAndScale>(reader.transformationToDensityLattice());
 153         }
 154         /*! \brief Set the periodic boundary condition via MdModuleNotifier.
 155          *
 156          * The pbc type is wrapped in PeriodicBoundaryConditionType to
 157          * allow the MdModuleNotifier to statically distinguish the callback
 158          * function type from other 'int' function callbacks.
 159          *
 160          * \param[in] pbc MdModuleNotification class that contains a variable
 161          *                that enumerates the periodic boundary condition.
 162          */
 163         void setPeriodicBoundaryConditionType(PeriodicBoundaryConditionType pbc)
 164         {
 165             pbcType_ = std::make_unique<int>(pbc.pbcType);
 166         }
 167
 168         //! Get the periodic boundary conditions
 169         int periodicBoundaryConditionType()
 170         {
 171             if (pbcType_ == nullptr)
 172             {
 173                 GMX_THROW(InternalError("Periodic boundary condition enum not set for density guided simulation."));
 174             }
 175             return *pbcType_;
 176         }
 177
 178     private:
 179         //! The reference density to fit to
 180         std::unique_ptr<MultiDimArray<std::vector<float>, dynamicExtents3D> > referenceDensity_;
 181         //! The coordinate transformation into the reference density
 182         std::unique_ptr<TranslateAndScale> transformationToDensityLattice_;
 183         //! The local atom set to act on
 184         std::unique_ptr<LocalAtomSet>      localAtomSet_;
 185         //! The type of periodic boundary conditions in the simulation
 186         std::unique_ptr<int>               pbcType_;
 187
 188         GMX_DISALLOW_COPY_AND_ASSIGN(DensityFittingSimulationParameterSetup);
 189 };
 190
 191 /*! \internal
 192  * \brief Density fitting
 193  *
 194  * Class that implements the density fitting forces and potential
 195  * \note the virial calculation is not yet implemented
 196  */
 197 class DensityFitting final : public IMDModule
 198 {
 199     public:
 200         /*! \brief Construct the density fitting module.
 201          *
 202          * \param[in] notifier allows the module to subscribe to notifications from MdModules.
 203          *
 204          * The density fitting code subscribes to these notifications:
 205          *   - setting atom group indices in the densityFittingOptions_ by
 206          *     taking a parmeter const IndexGroupsAndNames &
 207          *   - storing its internal parameters in a tpr file by writing to a
 208          *     key-value-tree during pre-processing by a function taking a
 209          *     KeyValueTreeObjectBuilder as parameter
 210          *   - reading its internal parameters from a key-value-tree during
 211          *     simulation setup by taking a const KeyValueTreeObject & parameter
 212          *   - constructing local atom sets in the simulation parameter setup
 213          *     by taking a LocalAtomSetManager * as parameter
 214          *   - the type of periodic boundary conditions that are used
 215          *     by taking a PeriodicBoundaryConditionTypeEnum as parameter
 216          */
 217         explicit DensityFitting(MdModulesNotifier *notifier)
 218         {
 219             // Callbacks for several kinds of MdModuleNotification are created
 220             // and subscribed, and will be dispatched correctly at run time
 221             // based on the type of the parameter required by the lambda.
 222
 223             // Setting atom group indices
 224             const auto setFitGroupIndicesFunction = [this](const IndexGroupsAndNames &indexGroupsAndNames) {
 225                     densityFittingOptions_.setFitGroupIndices(indexGroupsAndNames);
 226                 };
 227             notifier->notifier_.subscribe(setFitGroupIndicesFunction);
 228
 229             // Writing internal parameters during pre-processing
 230             const auto writeInternalParametersFunction = [this](KeyValueTreeObjectBuilder treeBuilder) {
 231                     densityFittingOptions_.writeInternalParametersToKvt(treeBuilder);
 232                 };
 233             notifier->notifier_.subscribe(writeInternalParametersFunction);
 234
 235             // Reading internal parameters during simulation setup
 236             const auto readInternalParametersFunction = [this](const KeyValueTreeObject &tree) {
 237                     densityFittingOptions_.readInternalParametersFromKvt(tree);
 238                 };
 239             notifier->notifier_.subscribe(readInternalParametersFunction);
 240
 241             // constructing local atom sets during simulation setup
 242             const auto setLocalAtomSetFunction = [this](LocalAtomSetManager *localAtomSetManager) {
 243                     this->constructLocalAtomSet(localAtomSetManager);
 244                 };
 245             notifier->notifier_.subscribe(setLocalAtomSetFunction);
 246
 247             // constructing local atom sets during simulation setup
 248             const auto setPeriodicBoundaryContionsFunction = [this](PeriodicBoundaryConditionType pbc) {
 249                     this->densityFittingSimulationParameters_.setPeriodicBoundaryConditionType(pbc);
 250                 };
 251             notifier->notifier_.subscribe(setPeriodicBoundaryContionsFunction);
 252         }
 253
 254         //! From IMDModule
 255         IMdpOptionProvider *mdpOptionProvider() override { return &densityFittingOptions_; }
 256
 257         //! Add this module to the force providers if active
 258         void initForceProviders(ForceProviders *forceProviders) override
 259         {
 260             if (densityFittingOptions_.active())
 261             {
 262                 const auto &parameters = densityFittingOptions_.buildParameters();
 263                 densityFittingSimulationParameters_.readReferenceDensityFromFile(densityFittingOptions_.referenceDensityFileName());
 264                 forceProvider_ = std::make_unique<DensityFittingForceProvider>(
 265                             parameters,
 266                             densityFittingSimulationParameters_.referenceDensity(),
 267                             densityFittingSimulationParameters_.transformationToDensityLattice(),
 268                             densityFittingSimulationParameters_.localAtomSet(),
 269                             densityFittingSimulationParameters_.periodicBoundaryConditionType());
 270                 forceProviders->addForceProvider(forceProvider_.get());
 271             }
 272         }
 273
 274         //! This MDModule provides its own output
 275         IMDOutputProvider *outputProvider() override { return &densityFittingOutputProvider_; }
 276
 277         /*! \brief Set up the local atom sets that are used by this module.
 278          *
 279          * \note When density fitting is set up with MdModuleNotification in
 280          *       the constructor, this function is called back.
 281          *
 282          * \param[in] localAtomSetManager the manager to add local atom sets.
 283          */
 284         void constructLocalAtomSet(LocalAtomSetManager * localAtomSetManager)
 285         {
 286             LocalAtomSet atomSet = localAtomSetManager->add(densityFittingOptions_.buildParameters().indices_);
 287             densityFittingSimulationParameters_.setLocalAtomSet(atomSet);
 288         }
 289
 290     private:
 291         //! The output provider
 292         DensityFittingOutputProvider                 densityFittingOutputProvider_;
 293         //! The options provided for density fitting
 294         DensityFittingOptions                        densityFittingOptions_;
 295         //! Object that evaluates the forces
 296         std::unique_ptr<DensityFittingForceProvider> forceProvider_;
 297         /*! \brief Parameters for density fitting that become available at
 298          * simulation setup time.
 299          */
 300         DensityFittingSimulationParameterSetup       densityFittingSimulationParameters_;
 301
 302         GMX_DISALLOW_COPY_AND_ASSIGN(DensityFitting);
 303 };
 304
 305 }   // namespace
 306
 307 std::unique_ptr<IMDModule> DensityFittingModuleInfo::create(MdModulesNotifier * notifier)
 308 {
 309     return std::make_unique<DensityFitting>(notifier);
 310 }
 311
 312 const std::string DensityFittingModuleInfo::name_ = "density-guided-simulation";
 313
 314 } // namespace gmx