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 #include <numeric>
  48
  49 #include "gromacs/domdec/localatomsetmanager.h"
  50 #include "gromacs/fileio/checkpoint.h"
  51 #include "gromacs/fileio/mrcdensitymap.h"
  52 #include "gromacs/math/multidimarray.h"
  53 #include "gromacs/mdlib/broadcaststructs.h"
  54 #include "gromacs/mdtypes/commrec.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(begin(referenceDensity_->asView()),
 167                                                       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] pbc MdModuleNotification class that contains a variable
 180      *                that enumerates the periodic boundary condition.
 181      */
 182     void setPeriodicBoundaryConditionType(PeriodicBoundaryConditionType pbc)
 183     {
 184         pbcType_ = std::make_unique<int>(pbc.pbcType);
 185     }
 186
 187     //! Get the periodic boundary conditions
 188     int periodicBoundaryConditionType()
 189     {
 190         if (pbcType_ == nullptr)
 191         {
 192             GMX_THROW(InternalError(
 193                     "Periodic boundary condition enum not set for density guided simulation."));
 194         }
 195         return *pbcType_;
 196     }
 197
 198     //! Set the simulation time step
 199     void setSimulationTimeStep(double timeStep) { simulationTimeStep_ = timeStep; }
 200
 201     //! Return the simulation time step
 202     double simulationTimeStep() { return simulationTimeStep_; }
 203
 204 private:
 205     //! The reference density to fit to
 206     std::unique_ptr<MultiDimArray<std::vector<float>, dynamicExtents3D>> referenceDensity_;
 207     //! The coordinate transformation into the reference density
 208     std::unique_ptr<TranslateAndScale> transformationToDensityLattice_;
 209     //! The local atom set to act on
 210     std::unique_ptr<LocalAtomSet> localAtomSet_;
 211     //! The type of periodic boundary conditions in the simulation
 212     std::unique_ptr<int> pbcType_;
 213     //! The simulation time step
 214     double simulationTimeStep_ = 1;
 215
 216     GMX_DISALLOW_COPY_AND_ASSIGN(DensityFittingSimulationParameterSetup);
 217 };
 218
 219 /*! \internal
 220  * \brief Density fitting
 221  *
 222  * Class that implements the density fitting forces and potential
 223  * \note the virial calculation is not yet implemented
 224  */
 225 class DensityFitting final : public IMDModule
 226 {
 227 public:
 228     /*! \brief Construct the density fitting module.
 229      *
 230      * \param[in] notifier allows the module to subscribe to notifications from MdModules.
 231      *
 232      * The density fitting code subscribes to these notifications:
 233      *   - setting atom group indices in the densityFittingOptions_ by
 234      *     taking a parmeter const IndexGroupsAndNames &
 235      *   - storing its internal parameters in a tpr file by writing to a
 236      *     key-value-tree during pre-processing by a function taking a
 237      *     KeyValueTreeObjectBuilder as parameter
 238      *   - reading its internal parameters from a key-value-tree during
 239      *     simulation setup by taking a const KeyValueTreeObject & parameter
 240      *   - constructing local atom sets in the simulation parameter setup
 241      *     by taking a LocalAtomSetManager * as parameter
 242      *   - the type of periodic boundary conditions that are used
 243      *     by taking a PeriodicBoundaryConditionType as parameter
 244      *   - the writing of checkpoint data
 245      *     by taking a MdModulesWriteCheckpointData as parameter
 246      *   - the reading of checkpoint data
 247      *     by taking a MdModulesCheckpointReadingDataOnMaster as parameter
 248      *   - the broadcasting of checkpoint data
 249      *     by taking MdModulesCheckpointReadingBroadcast as parameter
 250      */
 251     explicit DensityFitting(MdModulesNotifier* notifier)
 252     {
 253         // Callbacks for several kinds of MdModuleNotification are created
 254         // and subscribed, and will be dispatched correctly at run time
 255         // based on the type of the parameter required by the lambda.
 256
 257         // Setting atom group indices
 258         const auto setFitGroupIndicesFunction = [this](const IndexGroupsAndNames& indexGroupsAndNames) {
 259             densityFittingOptions_.setFitGroupIndices(indexGroupsAndNames);
 260         };
 261         notifier->notifier_.subscribe(setFitGroupIndicesFunction);
 262
 263         // Writing internal parameters during pre-processing
 264         const auto writeInternalParametersFunction = [this](KeyValueTreeObjectBuilder treeBuilder) {
 265             densityFittingOptions_.writeInternalParametersToKvt(treeBuilder);
 266         };
 267         notifier->notifier_.subscribe(writeInternalParametersFunction);
 268
 269         // Reading internal parameters during simulation setup
 270         const auto readInternalParametersFunction = [this](const KeyValueTreeObject& tree) {
 271             densityFittingOptions_.readInternalParametersFromKvt(tree);
 272         };
 273         notifier->notifier_.subscribe(readInternalParametersFunction);
 274
 275         // Checking for consistency with all .mdp options
 276         const auto checkEnergyCaluclationFrequencyFunction =
 277                 [this](EnergyCalculationFrequencyErrors* energyCalculationFrequencyErrors) {
 278                     densityFittingOptions_.checkEnergyCaluclationFrequency(energyCalculationFrequencyErrors);
 279                 };
 280         notifier->notifier_.subscribe(checkEnergyCaluclationFrequencyFunction);
 281
 282         // constructing local atom sets during simulation setup
 283         const auto setLocalAtomSetFunction = [this](LocalAtomSetManager* localAtomSetManager) {
 284             this->constructLocalAtomSet(localAtomSetManager);
 285         };
 286         notifier->notifier_.subscribe(setLocalAtomSetFunction);
 287
 288         // constructing local atom sets during simulation setup
 289         const auto setPeriodicBoundaryContionsFunction = [this](PeriodicBoundaryConditionType pbc) {
 290             this->densityFittingSimulationParameters_.setPeriodicBoundaryConditionType(pbc);
 291         };
 292         notifier->notifier_.subscribe(setPeriodicBoundaryContionsFunction);
 293
 294         // setting the simulation time step
 295         const auto setSimulationTimeStepFunction = [this](const SimulationTimeStep& simulationTimeStep) {
 296             this->densityFittingSimulationParameters_.setSimulationTimeStep(simulationTimeStep.delta_t);
 297         };
 298         notifier->notifier_.subscribe(setSimulationTimeStepFunction);
 299
 300         // adding output to energy file
 301         const auto requestEnergyOutput =
 302                 [this](MdModulesEnergyOutputToDensityFittingRequestChecker* energyOutputRequest) {
 303                     this->setEnergyOutputRequest(energyOutputRequest);
 304                 };
 305         notifier->notifier_.subscribe(requestEnergyOutput);
 306
 307         // writing checkpoint data
 308         const auto checkpointDataWriting = [this](MdModulesWriteCheckpointData checkpointData) {
 309             this->writeCheckpointData(checkpointData);
 310         };
 311         notifier->notifier_.subscribe(checkpointDataWriting);
 312
 313         // reading checkpoint data
 314         const auto checkpointDataReading = [this](MdModulesCheckpointReadingDataOnMaster checkpointData) {
 315             this->readCheckpointDataOnMaster(checkpointData);
 316         };
 317         notifier->notifier_.subscribe(checkpointDataReading);
 318
 319         // broadcasting checkpoint data
 320         const auto checkpointDataBroadcast = [this](MdModulesCheckpointReadingBroadcast checkpointData) {
 321             this->broadcastCheckpointData(checkpointData);
 322         };
 323         notifier->notifier_.subscribe(checkpointDataBroadcast);
 324     }
 325
 326     //! From IMDModule
 327     IMdpOptionProvider* mdpOptionProvider() override { return &densityFittingOptions_; }
 328
 329     //! Add this module to the force providers if active
 330     void initForceProviders(ForceProviders* forceProviders) override
 331     {
 332         if (densityFittingOptions_.active())
 333         {
 334             const auto& parameters = densityFittingOptions_.buildParameters();
 335             densityFittingSimulationParameters_.readReferenceDensityFromFile(
 336                     densityFittingOptions_.referenceDensityFileName());
 337             if (parameters.normalizeDensities_)
 338             {
 339                 densityFittingSimulationParameters_.normalizeReferenceDensity();
 340             }
 341             forceProvider_ = std::make_unique<DensityFittingForceProvider>(
 342                     parameters, densityFittingSimulationParameters_.referenceDensity(),
 343                     densityFittingSimulationParameters_.transformationToDensityLattice(),
 344                     densityFittingSimulationParameters_.localAtomSet(),
 345                     densityFittingSimulationParameters_.periodicBoundaryConditionType(),
 346                     densityFittingSimulationParameters_.simulationTimeStep(), densityFittingState_);
 347             forceProviders->addForceProvider(forceProvider_.get());
 348         }
 349     }
 350
 351     //! This MDModule provides its own output
 352     IMDOutputProvider* outputProvider() override { return &densityFittingOutputProvider_; }
 353
 354     /*! \brief Set up the local atom sets that are used by this module.
 355      *
 356      * \note When density fitting is set up with MdModuleNotification in
 357      *       the constructor, this function is called back.
 358      *
 359      * \param[in] localAtomSetManager the manager to add local atom sets.
 360      */
 361     void constructLocalAtomSet(LocalAtomSetManager* localAtomSetManager)
 362     {
 363         if (densityFittingOptions_.active())
 364         {
 365             LocalAtomSet atomSet =
 366                     localAtomSetManager->add(densityFittingOptions_.buildParameters().indices_);
 367             densityFittingSimulationParameters_.setLocalAtomSet(atomSet);
 368         }
 369     }
 370
 371     /*! \brief Request energy output to energy file during simulation.
 372      */
 373     void setEnergyOutputRequest(MdModulesEnergyOutputToDensityFittingRequestChecker* energyOutputRequest)
 374     {
 375         energyOutputRequest->energyOutputToDensityFitting_ = densityFittingOptions_.active();
 376     }
 377
 378     /*! \brief Write internal density fitting data to checkpoint file.
 379      * \param[in] checkpointWriting enables writing to the Key-Value-Tree
 380      *                              that is used for storing the checkpoint
 381      *                              information
 382      *
 383      * \note The provided state to checkpoint has to change if checkpointing
 384      *       is moved before the force provider call in the MD-loop.
 385      */
 386     void writeCheckpointData(MdModulesWriteCheckpointData checkpointWriting)
 387     {
 388         if (densityFittingOptions_.active())
 389         {
 390             const DensityFittingForceProviderState& state = forceProvider_->stateToCheckpoint();
 391             checkpointWriting.builder_.addValue<std::int64_t>(
 392                     DensityFittingModuleInfo::name_ + "-stepsSinceLastCalculation",
 393                     state.stepsSinceLastCalculation_);
 394             checkpointWriting.builder_.addValue<real>(
 395                     DensityFittingModuleInfo::name_ + "-adaptiveForceConstantScale",
 396                     state.adaptiveForceConstantScale_);
 397             KeyValueTreeObjectBuilder exponentialMovingAverageKvtEntry =
 398                     checkpointWriting.builder_.addObject(DensityFittingModuleInfo::name_
 399                                                          + "-exponentialMovingAverageState");
 400             exponentialMovingAverageStateAsKeyValueTree(exponentialMovingAverageKvtEntry,
 401                                                         state.exponentialMovingAverageState_);
 402         }
 403     }
 404
 405     /*! \brief Read the internal parameters from the checkpoint file on master
 406      * \param[in] checkpointReading holding the checkpoint information
 407      */
 408     void readCheckpointDataOnMaster(MdModulesCheckpointReadingDataOnMaster checkpointReading)
 409     {
 410         if (densityFittingOptions_.active())
 411         {
 412             if (checkpointReading.checkpointedData_.keyExists(DensityFittingModuleInfo::name_
 413                                                               + "-stepsSinceLastCalculation"))
 414             {
 415                 densityFittingState_.stepsSinceLastCalculation_ =
 416                         checkpointReading
 417                                 .checkpointedData_[DensityFittingModuleInfo::name_
 418                                                    + "-stepsSinceLastCalculation"]
 419                                 .cast<std::int64_t>();
 420             }
 421             if (checkpointReading.checkpointedData_.keyExists(DensityFittingModuleInfo::name_
 422                                                               + "-adaptiveForceConstantScale"))
 423             {
 424                 densityFittingState_.adaptiveForceConstantScale_ =
 425                         checkpointReading
 426                                 .checkpointedData_[DensityFittingModuleInfo::name_
 427                                                    + "-adaptiveForceConstantScale"]
 428                                 .cast<real>();
 429             }
 430             if (checkpointReading.checkpointedData_.keyExists(DensityFittingModuleInfo::name_
 431                                                               + "-exponentialMovingAverageState"))
 432             {
 433                 densityFittingState_.exponentialMovingAverageState_ = exponentialMovingAverageStateFromKeyValueTree(
 434                         checkpointReading
 435                                 .checkpointedData_[DensityFittingModuleInfo::name_ + "-exponentialMovingAverageState"]
 436                                 .asObject());
 437             }
 438         }
 439     }
 440
 441     /*! \brief Broadcast the internal parameters.
 442      * \param[in] checkpointBroadcast containing the communication record to
 443      *                                broadcast the checkpoint information
 444      */
 445     void broadcastCheckpointData(MdModulesCheckpointReadingBroadcast checkpointBroadcast)
 446     {
 447         if (densityFittingOptions_.active())
 448         {
 449             if (PAR(&(checkpointBroadcast.cr_)))
 450             {
 451                 block_bc(&(checkpointBroadcast.cr_), densityFittingState_.stepsSinceLastCalculation_);
 452                 block_bc(&(checkpointBroadcast.cr_), densityFittingState_.adaptiveForceConstantScale_);
 453                 block_bc(&(checkpointBroadcast.cr_), densityFittingState_.exponentialMovingAverageState_);
 454             }
 455         }
 456     }
 457
 458 private:
 459     //! The output provider
 460     DensityFittingOutputProvider densityFittingOutputProvider_;
 461     //! The options provided for density fitting
 462     DensityFittingOptions densityFittingOptions_;
 463     //! Object that evaluates the forces
 464     std::unique_ptr<DensityFittingForceProvider> forceProvider_;
 465     /*! \brief Parameters for density fitting that become available at
 466      * simulation setup time.
 467      */
 468     DensityFittingSimulationParameterSetup densityFittingSimulationParameters_;
 469     //! The internal parameters of density fitting force provider
 470     DensityFittingForceProviderState densityFittingState_;
 471
 472     GMX_DISALLOW_COPY_AND_ASSIGN(DensityFitting);
 473 };
 474
 475 } // namespace
 476
 477 std::unique_ptr<IMDModule> DensityFittingModuleInfo::create(MdModulesNotifier* notifier)
 478 {
 479     return std::make_unique<DensityFitting>(notifier);
 480 }
 481
 482 const std::string DensityFittingModuleInfo::name_ = "density-guided-simulation";
 483
 484 } // namespace gmx