src/gromacs/applied_forces/densityfitting/densityfitting.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020, 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 #include <numeric>
  48
  49 #include "gromacs/domdec/localatomset.h"
  50 #include "gromacs/domdec/localatomsetmanager.h"
  51 #include "gromacs/fileio/checkpoint.h"
  52 #include "gromacs/fileio/mrcdensitymap.h"
  53 #include "gromacs/math/coordinatetransformation.h"
  54 #include "gromacs/math/multidimarray.h"
  55 #include "gromacs/mdtypes/imdmodule.h"
  56 #include "gromacs/utility/classhelpers.h"
  57 #include "gromacs/utility/exceptions.h"
  58 #include "gromacs/utility/keyvaluetreebuilder.h"
  59 #include "gromacs/utility/mdmodulenotification.h"
  60
  61 #include "densityfittingforceprovider.h"
  62 #include "densityfittingoptions.h"
  63 #include "densityfittingoutputprovider.h"
  64
  65 namespace gmx
  66 {
  67
  68 class IMdpOptionProvider;
  69 class DensityFittingForceProvider;
  70
  71 namespace
  72 {
  73
  74 /*! \internal
  75  * \brief Collect density fitting parameters only available during simulation setup.
  76  *
  77  * \todo Implement builder pattern that will not use unqiue_ptr to check if
  78  *       parameters have been set or not.
  79  *
  80  * To build the density fitting force provider during simulation setup,
  81  * the DensityFitting class needs access to parameters that become available
  82  * only during simulation setup.
  83  *
  84  * This class collects these parameters via MdModuleNotifications in the
  85  * simulation setup phase and provides a check if all necessary parameters have
  86  * been provided.
  87  */
  88 class DensityFittingSimulationParameterSetup
  89 {
  90 public:
  91     DensityFittingSimulationParameterSetup() = default;
  92
  93     /*! \brief Set the local atom set for the density fitting.
  94      * \param[in] localAtomSet of atoms to be fitted
  95      */
  96     void setLocalAtomSet(const LocalAtomSet& localAtomSet)
  97     {
  98         localAtomSet_ = std::make_unique<LocalAtomSet>(localAtomSet);
  99     }
 100
 101     /*! \brief Return local atom set for density fitting.
 102      * \throws InternalError if local atom set is not set
 103      * \returns local atom set for density fitting
 104      */
 105     const LocalAtomSet& localAtomSet() const
 106     {
 107         if (localAtomSet_ == nullptr)
 108         {
 109             GMX_THROW(
 110                     InternalError("Local atom set is not set for density "
 111                                   "guided simulation."));
 112         }
 113         return *localAtomSet_;
 114     }
 115
 116     /*! \brief Return transformation into density lattice.
 117      * \throws InternalError if transformation into density lattice is not set
 118      * \returns transformation into density lattice
 119      */
 120     const TranslateAndScale& transformationToDensityLattice() const
 121     {
 122         if (transformationToDensityLattice_ == nullptr)
 123         {
 124             GMX_THROW(InternalError(
 125                     "Transformation to reference density not set for density guided simulation."));
 126         }
 127         return *transformationToDensityLattice_;
 128     }
 129     /*! \brief Return reference density
 130      * \throws InternalError if reference density is not set
 131      * \returns the reference density
 132      */
 133     basic_mdspan<const float, dynamicExtents3D> referenceDensity() const
 134     {
 135         if (referenceDensity_ == nullptr)
 136         {
 137             GMX_THROW(InternalError("Reference density not set for density guided simulation."));
 138         }
 139         return referenceDensity_->asConstView();
 140     }
 141
 142     /*! \brief Reads the reference density from file.
 143      *
 144      * Reads and check file, then set and communicate the internal
 145      * parameters related to the reference density with the file data.
 146      *
 147      * \throws FileIOError if reading from file was not successful
 148      */
 149     void readReferenceDensityFromFile(const std::string& referenceDensityFileName)
 150     {
 151         MrcDensityMapOfFloatFromFileReader reader(referenceDensityFileName);
 152         referenceDensity_ = std::make_unique<MultiDimArray<std::vector<float>, dynamicExtents3D>>(
 153                 reader.densityDataCopy());
 154         transformationToDensityLattice_ =
 155                 std::make_unique<TranslateAndScale>(reader.transformationToDensityLattice());
 156     }
 157
 158     //! Normalize the reference density so that the sum over all voxels is unity
 159     void normalizeReferenceDensity()
 160     {
 161         if (referenceDensity_ == nullptr)
 162         {
 163             GMX_THROW(InternalError("Need to set reference density before normalizing it."));
 164         }
 165
 166         const real sumOfDensityData = std::accumulate(
 167                 begin(referenceDensity_->asView()), end(referenceDensity_->asView()), 0.);
 168         for (float& referenceDensityVoxel : referenceDensity_->asView())
 169         {
 170             referenceDensityVoxel /= sumOfDensityData;
 171         }
 172     }
 173     /*! \brief Set the periodic boundary condition via MdModuleNotifier.
 174      *
 175      * The pbc type is wrapped in PeriodicBoundaryConditionType to
 176      * allow the MdModuleNotifier to statically distinguish the callback
 177      * function type from other 'int' function callbacks.
 178      *
 179      * \param[in] pbcType enumerates the periodic boundary condition.
 180      */
 181     void setPeriodicBoundaryConditionType(const PbcType& pbcType)
 182     {
 183         pbcType_ = std::make_unique<PbcType>(pbcType);
 184     }
 185
 186     //! Get the periodic boundary conditions
 187     PbcType periodicBoundaryConditionType()
 188     {
 189         if (pbcType_ == nullptr)
 190         {
 191             GMX_THROW(InternalError(
 192                     "Periodic boundary condition enum not set for density guided simulation."));
 193         }
 194         return *pbcType_;
 195     }
 196
 197     //! Set the simulation time step
 198     void setSimulationTimeStep(double timeStep) { simulationTimeStep_ = timeStep; }
 199
 200     //! Return the simulation time step
 201     double simulationTimeStep() { return simulationTimeStep_; }
 202
 203 private:
 204     //! The reference density to fit to
 205     std::unique_ptr<MultiDimArray<std::vector<float>, dynamicExtents3D>> referenceDensity_;
 206     //! The coordinate transformation into the reference density
 207     std::unique_ptr<TranslateAndScale> transformationToDensityLattice_;
 208     //! The local atom set to act on
 209     std::unique_ptr<LocalAtomSet> localAtomSet_;
 210     //! The type of periodic boundary conditions in the simulation
 211     std::unique_ptr<PbcType> pbcType_;
 212     //! The simulation time step
 213     double simulationTimeStep_ = 1;
 214
 215     GMX_DISALLOW_COPY_AND_ASSIGN(DensityFittingSimulationParameterSetup);
 216 };
 217
 218 /*! \internal
 219  * \brief Density fitting
 220  *
 221  * Class that implements the density fitting forces and potential
 222  * \note the virial calculation is not yet implemented
 223  */
 224 class DensityFitting final : public IMDModule
 225 {
 226
 227 public:
 228     //! Construct the density fitting module.
 229     explicit DensityFitting() = default;
 230
 231     /*! \brief Request to be notified during pre-processing.
 232      *
 233      * \param[in] notifier allows the module to subscribe to notifications from MdModules.
 234      *
 235      * The density fitting code subscribes to these notifications:
 236      *   - setting atom group indices in the densityFittingOptions_ from an
 237      *     index group string by taking a parmeter const IndexGroupsAndNames &
 238      *   - storing its internal parameters in a tpr file by writing to a
 239      *     key-value-tree during pre-processing by a function taking a
 240      *     KeyValueTreeObjectBuilder as parameter
 241      */
 242     void subscribeToPreProcessingNotifications(MdModulesNotifier* notifier) override
 243     {
 244         // Callbacks for several kinds of MdModuleNotification are created
 245         // and subscribed, and will be dispatched correctly at run time
 246         // based on the type of the parameter required by the lambda.
 247
 248         if (!densityFittingOptions_.active())
 249         {
 250             return;
 251         }
 252
 253         // Setting the atom group indices from index group string
 254         const auto setFitGroupIndicesFunction = [this](const IndexGroupsAndNames& indexGroupsAndNames) {
 255             densityFittingOptions_.setFitGroupIndices(indexGroupsAndNames);
 256         };
 257         notifier->preProcessingNotifications_.subscribe(setFitGroupIndicesFunction);
 258
 259         // Writing internal parameters during pre-processing
 260         const auto writeInternalParametersFunction = [this](KeyValueTreeObjectBuilder treeBuilder) {
 261             densityFittingOptions_.writeInternalParametersToKvt(treeBuilder);
 262         };
 263         notifier->preProcessingNotifications_.subscribe(writeInternalParametersFunction);
 264
 265         // Checking for consistency with all .mdp options
 266         const auto checkEnergyCaluclationFrequencyFunction =
 267                 [this](EnergyCalculationFrequencyErrors* energyCalculationFrequencyErrors) {
 268                     densityFittingOptions_.checkEnergyCaluclationFrequency(energyCalculationFrequencyErrors);
 269                 };
 270         notifier->preProcessingNotifications_.subscribe(checkEnergyCaluclationFrequencyFunction);
 271     }
 272
 273     /*! \brief Request to be notified.
 274      * The density fitting code subscribes to these notifications:
 275      *   - reading its internal parameters from a key-value-tree during
 276      *     simulation setup by taking a const KeyValueTreeObject & parameter
 277      *   - constructing local atom sets in the simulation parameter setup
 278      *     by taking a LocalAtomSetManager * as parameter
 279      *   - the type of periodic boundary conditions that are used
 280      *     by taking a PeriodicBoundaryConditionType as parameter
 281      *   - the writing of checkpoint data
 282      *     by taking a MdModulesWriteCheckpointData as parameter
 283      *   - the reading of checkpoint data
 284      *     by taking a MdModulesCheckpointReadingDataOnMaster as parameter
 285      *   - the broadcasting of checkpoint data
 286      *     by taking MdModulesCheckpointReadingBroadcast as parameter
 287      */
 288     void subscribeToSimulationSetupNotifications(MdModulesNotifier* notifier) override
 289     {
 290         if (!densityFittingOptions_.active())
 291         {
 292             return;
 293         }
 294
 295         // Reading internal parameters during simulation setup
 296         const auto readInternalParametersFunction = [this](const KeyValueTreeObject& tree) {
 297             densityFittingOptions_.readInternalParametersFromKvt(tree);
 298         };
 299         notifier->simulationSetupNotifications_.subscribe(readInternalParametersFunction);
 300
 301         // constructing local atom sets during simulation setup
 302         const auto setLocalAtomSetFunction = [this](LocalAtomSetManager* localAtomSetManager) {
 303             this->constructLocalAtomSet(localAtomSetManager);
 304         };
 305         notifier->simulationSetupNotifications_.subscribe(setLocalAtomSetFunction);
 306
 307         // constructing local atom sets during simulation setup
 308         const auto setPeriodicBoundaryContionsFunction = [this](const PbcType& pbc) {
 309             this->densityFittingSimulationParameters_.setPeriodicBoundaryConditionType(pbc);
 310         };
 311         notifier->simulationSetupNotifications_.subscribe(setPeriodicBoundaryContionsFunction);
 312
 313         // setting the simulation time step
 314         const auto setSimulationTimeStepFunction = [this](const SimulationTimeStep& simulationTimeStep) {
 315             this->densityFittingSimulationParameters_.setSimulationTimeStep(simulationTimeStep.delta_t);
 316         };
 317         notifier->simulationSetupNotifications_.subscribe(setSimulationTimeStepFunction);
 318
 319         // adding output to energy file
 320         const auto requestEnergyOutput =
 321                 [this](MdModulesEnergyOutputToDensityFittingRequestChecker* energyOutputRequest) {
 322                     this->setEnergyOutputRequest(energyOutputRequest);
 323                 };
 324         notifier->simulationSetupNotifications_.subscribe(requestEnergyOutput);
 325
 326         // writing checkpoint data
 327         const auto checkpointDataWriting = [this](MdModulesWriteCheckpointData checkpointData) {
 328             forceProvider_->writeCheckpointData(checkpointData, DensityFittingModuleInfo::name_);
 329         };
 330         notifier->checkpointingNotifications_.subscribe(checkpointDataWriting);
 331
 332         // reading checkpoint data
 333         const auto checkpointDataReading = [this](MdModulesCheckpointReadingDataOnMaster checkpointData) {
 334             densityFittingState_.readState(checkpointData.checkpointedData_,
 335                                            DensityFittingModuleInfo::name_);
 336         };
 337         notifier->checkpointingNotifications_.subscribe(checkpointDataReading);
 338
 339         // broadcasting checkpoint data
 340         const auto checkpointDataBroadcast = [this](MdModulesCheckpointReadingBroadcast checkpointData) {
 341             densityFittingState_.broadcastState(checkpointData.communicator_, checkpointData.isParallelRun_);
 342         };
 343         notifier->checkpointingNotifications_.subscribe(checkpointDataBroadcast);
 344     }
 345
 346     //! From IMDModule
 347     IMdpOptionProvider* mdpOptionProvider() override { return &densityFittingOptions_; }
 348
 349     //! Add this module to the force providers if active
 350     void initForceProviders(ForceProviders* forceProviders) override
 351     {
 352         if (densityFittingOptions_.active())
 353         {
 354             const auto& parameters = densityFittingOptions_.buildParameters();
 355             densityFittingSimulationParameters_.readReferenceDensityFromFile(
 356                     densityFittingOptions_.referenceDensityFileName());
 357             if (parameters.normalizeDensities_)
 358             {
 359                 densityFittingSimulationParameters_.normalizeReferenceDensity();
 360             }
 361             forceProvider_ = std::make_unique<DensityFittingForceProvider>(
 362                     parameters,
 363                     densityFittingSimulationParameters_.referenceDensity(),
 364                     densityFittingSimulationParameters_.transformationToDensityLattice(),
 365                     densityFittingSimulationParameters_.localAtomSet(),
 366                     densityFittingSimulationParameters_.periodicBoundaryConditionType(),
 367                     densityFittingSimulationParameters_.simulationTimeStep(),
 368                     densityFittingState_);
 369             forceProviders->addForceProvider(forceProvider_.get());
 370         }
 371     }
 372
 373     //! This MDModule provides its own output
 374     IMDOutputProvider* outputProvider() override { return &densityFittingOutputProvider_; }
 375
 376     /*! \brief Set up the local atom sets that are used by this module.
 377      *
 378      * \note When density fitting is set up with MdModuleNotification in
 379      *       the constructor, this function is called back.
 380      *
 381      * \param[in] localAtomSetManager the manager to add local atom sets.
 382      */
 383     void constructLocalAtomSet(LocalAtomSetManager* localAtomSetManager)
 384     {
 385         LocalAtomSet atomSet = localAtomSetManager->add(densityFittingOptions_.buildParameters().indices_);
 386         densityFittingSimulationParameters_.setLocalAtomSet(atomSet);
 387     }
 388
 389     /*! \brief Request energy output to energy file during simulation.
 390      */
 391     void setEnergyOutputRequest(MdModulesEnergyOutputToDensityFittingRequestChecker* energyOutputRequest)
 392     {
 393         energyOutputRequest->energyOutputToDensityFitting_ = densityFittingOptions_.active();
 394     }
 395
 396 private:
 397     //! The output provider
 398     DensityFittingOutputProvider densityFittingOutputProvider_;
 399     //! The options provided for density fitting
 400     DensityFittingOptions densityFittingOptions_;
 401     //! Object that evaluates the forces
 402     std::unique_ptr<DensityFittingForceProvider> forceProvider_;
 403     /*! \brief Parameters for density fitting that become available at
 404      * simulation setup time.
 405      */
 406     DensityFittingSimulationParameterSetup densityFittingSimulationParameters_;
 407     //! The internal parameters of density fitting force provider
 408     DensityFittingForceProviderState densityFittingState_;
 409
 410     GMX_DISALLOW_COPY_AND_ASSIGN(DensityFitting);
 411 };
 412
 413 } // namespace
 414
 415 std::unique_ptr<IMDModule> DensityFittingModuleInfo::create()
 416 {
 417     return std::make_unique<DensityFitting>();
 418 }
 419
 420 const std::string DensityFittingModuleInfo::name_ = "density-guided-simulation";
 421
 422 } // namespace gmx