src/gromacs/applied_forces/qmmm/qmmm.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2021, 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  * Implements QMMM class that implements IMDModule interface
  38  *
  39  * \author Dmitry Morozov <dmitry.morozov@jyu.fi>
  40  * \author Christian Blau <blau@kth.se>
  41  * \ingroup module_applied_forces
  42  */
  43 #include "gmxpre.h"
  44
  45 #include "qmmm.h"
  46
  47 #include <memory>
  48 #include <numeric>
  49
  50 #include "gromacs/domdec/localatomsetmanager.h"
  51 #include "gromacs/ewald/pme.h"
  52 #include "gromacs/fileio/checkpoint.h"
  53 #include "gromacs/math/multidimarray.h"
  54 #include "gromacs/mdlib/broadcaststructs.h"
  55 #include "gromacs/mdtypes/commrec.h"
  56 #include "gromacs/mdtypes/imdmodule.h"
  57 #include "gromacs/mdtypes/imdoutputprovider.h"
  58 #include "gromacs/utility/classhelpers.h"
  59 #include "gromacs/utility/exceptions.h"
  60 #include "gromacs/utility/keyvaluetreebuilder.h"
  61 #include "gromacs/utility/logger.h"
  62 #include "gromacs/utility/mdmodulesnotifiers.h"
  63
  64 #include "qmmmforceprovider.h"
  65 #include "qmmmoptions.h"
  66
  67 namespace gmx
  68 {
  69
  70 namespace
  71 {
  72
  73 /*! \internal
  74  * \brief Helper class that holds simulation data and
  75  * callback functions for simulation setup time notifications
  76  */
  77 class QMMMSimulationParameterSetup
  78 {
  79 public:
  80     QMMMSimulationParameterSetup() = default;
  81
  82     /*! \brief Set the local atom set for the QM atoms.
  83      * \param[in] localAtomSet of QM atoms
  84      */
  85     void setLocalQMAtomSet(const LocalAtomSet& localAtomSet)
  86     {
  87         localQMAtomSet_ = std::make_unique<LocalAtomSet>(localAtomSet);
  88     }
  89
  90     /*! \brief Set the local atom set for the MM atoms.
  91      * \param[in] localAtomSet of MM atoms
  92      */
  93     void setLocalMMAtomSet(const LocalAtomSet& localAtomSet)
  94     {
  95         localMMAtomSet_ = std::make_unique<LocalAtomSet>(localAtomSet);
  96     }
  97
  98     /*! \brief Return local atom set for QM atoms.
  99      * \throws InternalError if local atom set is not set
 100      * \returns local atom set for QM atoms
 101      */
 102     const LocalAtomSet& localQMAtomSet() const
 103     {
 104         if (localQMAtomSet_ == nullptr)
 105         {
 106             GMX_THROW(InternalError("Local atom set is not set for QM atoms."));
 107         }
 108         return *localQMAtomSet_;
 109     }
 110
 111     /*! \brief Return local atom set for MM atoms.
 112      * \throws InternalError if local atom set is not set
 113      * \returns local atom set for MM atoms
 114      */
 115     const LocalAtomSet& localMMAtomSet() const
 116     {
 117         if (localMMAtomSet_ == nullptr)
 118         {
 119             GMX_THROW(InternalError("Local atom set is not set for MM atoms."));
 120         }
 121         return *localMMAtomSet_;
 122     }
 123
 124     /*! \brief Set the periodic boundary condition via MdModuleNotifier.
 125      *
 126      * The pbc type is wrapped in PeriodicBoundaryConditionType to
 127      * allow the MdModuleNotifier to statically distinguish the callback
 128      * function type from other 'int' function callbacks.
 129      *
 130      * \param[in] pbc MdModuleNotification class that contains a variable
 131      *                that enumerates the periodic boundary condition.
 132      */
 133     void setPeriodicBoundaryConditionType(const PbcType& pbc)
 134     {
 135         pbcType_ = std::make_unique<PbcType>(pbc);
 136     }
 137
 138     //! Get the periodic boundary conditions
 139     PbcType periodicBoundaryConditionType()
 140     {
 141         if (pbcType_ == nullptr)
 142         {
 143             GMX_THROW(
 144                     InternalError("Periodic boundary condition enum not set for QMMM simulation."));
 145         }
 146         return *pbcType_;
 147     }
 148
 149     /*! \brief Set the logger for QMMM during mdrun
 150      * \param[in] logger Logger instance to be used for output
 151      */
 152     void setLogger(const MDLogger& logger) { logger_ = &logger; }
 153
 154     //! Get the logger instance
 155     const MDLogger& logger() const
 156     {
 157         if (logger_ == nullptr)
 158         {
 159             GMX_THROW(InternalError("Logger not set for QMMM simulation."));
 160         }
 161         return *logger_;
 162     }
 163
 164 private:
 165     //! The local QM atom set to act on
 166     std::unique_ptr<LocalAtomSet> localQMAtomSet_;
 167     //! The local MM atom set to act on
 168     std::unique_ptr<LocalAtomSet> localMMAtomSet_;
 169     //! The type of periodic boundary conditions in the simulation
 170     std::unique_ptr<PbcType> pbcType_;
 171     //! MDLogger for notifications during mdrun
 172     const MDLogger* logger_ = nullptr;
 173
 174     GMX_DISALLOW_COPY_AND_ASSIGN(QMMMSimulationParameterSetup);
 175 };
 176
 177
 178 /*! \internal
 179  * \brief Handle file output for QMMM simulations.
 180  * empty implementation as QMMM does not use that
 181  */
 182 class QMMMOutputProvider final : public IMDOutputProvider
 183 {
 184 public:
 185     //! Initialize output
 186     void initOutput(FILE* /*fplog*/,
 187                     int /*nfile*/,
 188                     const t_filenm /*fnm*/[],
 189                     bool /*bAppendFiles*/,
 190                     const gmx_output_env_t* /*oenv*/) override
 191     {
 192     }
 193     //! Finalizes output from a simulation run.
 194     void finishOutput() override {}
 195 };
 196
 197
 198 /*! \internal
 199  * \brief QMMM module
 200  *
 201  * Class that implements the QMMM forces and potential
 202  * \note the virial calculation is not yet implemented
 203  */
 204 class QMMM final : public IMDModule
 205 {
 206 public:
 207     //! \brief Construct the QMMM module.
 208     explicit QMMM() = default;
 209
 210     // Now callbacks for several kinds of MdModuleNotification are created
 211     // and subscribed, and will be dispatched correctly at run time
 212     // based on the type of the parameter required by the lambda.
 213
 214     /*! \brief Requests to be notified during pre-processing.
 215      *
 216      * \param[in] notifier allows the module to subscribe to notifications from MdModules.
 217      *
 218      * The QMMM code subscribes to these notifications:
 219      *   - setting atom group indices in the qmmmOptions_ from an
 220      *     index group string by taking a parmeter const IndexGroupsAndNames &
 221      *   - storing its internal parameters in a tpr file by writing to a
 222      *     key-value-tree during pre-processing by a function taking a
 223      *     KeyValueTreeObjectBuilder as parameter
 224      *   - Modify topology according to QMMM rules using gmx_mtop_t notification
 225      *     and utilizing QMMMTopologyPreprocessor class
 226      *   - Access MDLogger for notifications output
 227      *   - Access warninp for for grompp warnings output
 228      *   - Coordinates, PBC and box for CP2K input generation
 229      *   - QM Input file provided with -qmi option of grompp
 230      */
 231     void subscribeToPreProcessingNotifications(MDModulesNotifiers* notifier) override
 232     {
 233         if (!qmmmOptions_.active())
 234         {
 235             return;
 236         }
 237
 238         // Writing internal parameters during pre-processing
 239         const auto writeInternalParametersFunction = [this](KeyValueTreeObjectBuilder treeBuilder) {
 240             qmmmOptions_.writeInternalParametersToKvt(treeBuilder);
 241         };
 242         notifier->preProcessingNotifier_.subscribe(writeInternalParametersFunction);
 243
 244         // Setting atom group indices
 245         const auto setQMMMGroupIndicesFunction = [this](const IndexGroupsAndNames& indexGroupsAndNames) {
 246             qmmmOptions_.setQMMMGroupIndices(indexGroupsAndNames);
 247         };
 248         notifier->preProcessingNotifier_.subscribe(setQMMMGroupIndicesFunction);
 249
 250         // Set Logger during pre-processing
 251         const auto setLoggerFunction = [this](const MDLogger& logger) {
 252             qmmmOptions_.setLogger(logger);
 253         };
 254         notifier->preProcessingNotifier_.subscribe(setLoggerFunction);
 255
 256         // Set warning output during pre-processing
 257         const auto setWarninpFunction = [this](warninp* wi) { qmmmOptions_.setWarninp(wi); };
 258         notifier->preProcessingNotifier_.subscribe(setWarninpFunction);
 259
 260         // Notification of the Coordinates, box and pbc during pre-processing
 261         const auto processCoordinatesFunction = [this](const CoordinatesAndBoxPreprocessed& coord) {
 262             qmmmOptions_.processCoordinates(coord);
 263         };
 264         notifier->preProcessingNotifier_.subscribe(processCoordinatesFunction);
 265
 266         // Modification of the topology during pre-processing
 267         const auto modifyQMMMTopologyFunction = [this](gmx_mtop_t* mtop) {
 268             qmmmOptions_.modifyQMMMTopology(mtop);
 269         };
 270         notifier->preProcessingNotifier_.subscribe(modifyQMMMTopologyFunction);
 271
 272         // Notification of the QM input file provided via -qmi option of grompp
 273         const auto setQMExternalInputFileNameFunction = [this](const QMInputFileName& qmInputFileName) {
 274             qmmmOptions_.setQMExternalInputFile(qmInputFileName);
 275         };
 276         notifier->preProcessingNotifier_.subscribe(setQMExternalInputFileNameFunction);
 277     }
 278
 279     /*! \brief Requests to be notified during simulation setup.
 280      * The QMMM code subscribes to these notifications:
 281      *   - reading its internal parameters from a key-value-tree during
 282      *     simulation setup by taking a const KeyValueTreeObject & parameter
 283      *   - *.tpr filename for CP2K input generation
 284      *   - constructing local atom sets in the simulation parameter setup
 285      *     by taking a LocalAtomSetManager * as parameter
 286      *   - the type of periodic boundary conditions that are used
 287      *     by taking a PeriodicBoundaryConditionType as parameter
 288      *   - Access MDLogger for notifications output
 289      *   - Disable PME-only ranks for QMMM runs
 290      *   - Request QM energy output to md.log
 291      */
 292     void subscribeToSimulationSetupNotifications(MDModulesNotifiers* notifier) override
 293     {
 294         if (!qmmmOptions_.active())
 295         {
 296             return;
 297         }
 298
 299         // Reading internal parameters during simulation setup
 300         const auto readInternalParametersFunction = [this](const KeyValueTreeObject& tree) {
 301             qmmmOptions_.readInternalParametersFromKvt(tree);
 302         };
 303         notifier->simulationSetupNotifier_.subscribe(readInternalParametersFunction);
 304
 305         // Process tpr filename
 306         const auto setTprFileNameFunction = [this](const MdRunInputFilename& tprName) {
 307             qmmmOptions_.processTprFilename(tprName);
 308         };
 309         notifier->simulationSetupNotifier_.subscribe(setTprFileNameFunction);
 310
 311         // constructing local atom sets during simulation setup
 312         const auto setLocalAtomSetFunction = [this](LocalAtomSetManager* localAtomSetManager) {
 313             LocalAtomSet atomSet1 = localAtomSetManager->add(qmmmOptions_.parameters().qmIndices_);
 314             this->qmmmSimulationParameters_.setLocalQMAtomSet(atomSet1);
 315             LocalAtomSet atomSet2 = localAtomSetManager->add(qmmmOptions_.parameters().mmIndices_);
 316             this->qmmmSimulationParameters_.setLocalMMAtomSet(atomSet2);
 317         };
 318         notifier->simulationSetupNotifier_.subscribe(setLocalAtomSetFunction);
 319
 320         // Reading PBC parameters during simulation setup
 321         const auto setPeriodicBoundaryContionsFunction = [this](const PbcType& pbc) {
 322             this->qmmmSimulationParameters_.setPeriodicBoundaryConditionType(pbc);
 323         };
 324         notifier->simulationSetupNotifier_.subscribe(setPeriodicBoundaryContionsFunction);
 325
 326         // Saving MDLogger during simulation setup
 327         const auto setLoggerFunction = [this](const MDLogger& logger) {
 328             this->qmmmSimulationParameters_.setLogger(logger);
 329         };
 330         notifier->simulationSetupNotifier_.subscribe(setLoggerFunction);
 331
 332         // Adding output to energy file
 333         const auto requestEnergyOutput = [](MDModulesEnergyOutputToQMMMRequestChecker* energyOutputRequest) {
 334             energyOutputRequest->energyOutputToQMMM_ = true;
 335         };
 336         notifier->simulationSetupNotifier_.subscribe(requestEnergyOutput);
 337
 338         // Request to disable PME-only ranks, which are not compatible with CP2K
 339         const auto requestPmeRanks = [](SeparatePmeRanksPermitted* pmeRanksPermitted) {
 340             pmeRanksPermitted->disablePmeRanks(
 341                     "Separate PME-only ranks are not compatible with QMMM MdModule");
 342         };
 343         notifier->simulationSetupNotifier_.subscribe(requestPmeRanks);
 344     }
 345
 346     //! From IMDModule
 347     IMdpOptionProvider* mdpOptionProvider() override { return &qmmmOptions_; }
 348
 349     //! Add this module to the force providers if active
 350     void initForceProviders(ForceProviders* forceProviders) override
 351     {
 352         if (!qmmmOptions_.active())
 353         {
 354             return;
 355         }
 356
 357         const auto& parameters = qmmmOptions_.parameters();
 358         forceProvider_         = std::make_unique<QMMMForceProvider>(
 359                 parameters,
 360                 qmmmSimulationParameters_.localQMAtomSet(),
 361                 qmmmSimulationParameters_.localMMAtomSet(),
 362                 qmmmSimulationParameters_.periodicBoundaryConditionType(),
 363                 qmmmSimulationParameters_.logger());
 364         forceProviders->addForceProvider(forceProvider_.get());
 365     }
 366
 367     //! QMMM Module should not use OutputProvider as it will be removed in the furture
 368     IMDOutputProvider* outputProvider() override { return &qmmmOutputProvider_; }
 369
 370 private:
 371     //! The output provider
 372     QMMMOutputProvider qmmmOutputProvider_;
 373     //! The options provided for QMMM
 374     QMMMOptions qmmmOptions_;
 375     //! Object that evaluates the forces
 376     std::unique_ptr<QMMMForceProvider> forceProvider_;
 377     /*! \brief Parameters for QMMM that become available at
 378      * simulation setup time.
 379      */
 380     QMMMSimulationParameterSetup qmmmSimulationParameters_;
 381
 382     GMX_DISALLOW_COPY_AND_ASSIGN(QMMM);
 383 };
 384
 385 } // namespace
 386
 387 std::unique_ptr<IMDModule> QMMMModuleInfo::create()
 388 {
 389     return std::make_unique<QMMM>();
 390 }
 391
 392 const std::string QMMMModuleInfo::name_ = c_qmmmCP2KModuleName;
 393
 394 } // namespace gmx