src/gromacs/modularsimulator/firstorderpressurecoupling.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 Defines the element performing first-order pressure coupling for the modular simulator
  37  *
  38  * \author Pascal Merz <pascal.merz@me.com>
  39  * \ingroup module_modularsimulator
  40  */
  41
  42 #include "gmxpre.h"
  43
  44 #include "firstorderpressurecoupling.h"
  45
  46 #include "gromacs/domdec/domdec_network.h"
  47 #include "gromacs/mdlib/coupling.h"
  48 #include "gromacs/mdlib/mdatoms.h"
  49 #include "gromacs/mdlib/stat.h"
  50 #include "gromacs/mdtypes/commrec.h"
  51 #include "gromacs/mdtypes/inputrec.h"
  52 #include "gromacs/mdtypes/md_enums.h"
  53 #include "gromacs/mdtypes/mdatom.h"
  54 #include "gromacs/pbcutil/boxutilities.h"
  55
  56 #include "energydata.h"
  57 #include "statepropagatordata.h"
  58 #include "simulatoralgorithm.h"
  59
  60 namespace gmx
  61 {
  62
  63 template<PressureCoupling pressureCouplingType>
  64 void FirstOrderPressureCoupling::calculateScalingMatrix(Step step)
  65 {
  66     const auto* pressure         = energyData_->pressure(step);
  67     const auto* forceVirial      = energyData_->forceVirial(step);
  68     const auto* constraintVirial = energyData_->constraintVirial(step);
  69     const auto* box              = statePropagatorData_->constBox();
  70
  71     previousStepConservedEnergyContribution_ = conservedEnergyContribution_;
  72     pressureCouplingCalculateScalingMatrix<pressureCouplingType>(fplog_,
  73                                                                  step,
  74                                                                  inputrec_,
  75                                                                  couplingTimeStep_,
  76                                                                  pressure,
  77                                                                  box,
  78                                                                  forceVirial,
  79                                                                  constraintVirial,
  80                                                                  boxScalingMatrix_,
  81                                                                  &conservedEnergyContribution_);
  82     conservedEnergyContributionStep_ = step;
  83 }
  84
  85 template<PressureCoupling pressureCouplingType>
  86 void FirstOrderPressureCoupling::scaleBoxAndCoordinates()
  87 {
  88     auto*          box       = statePropagatorData_->box();
  89     auto           positions = statePropagatorData_->positionsView().unpaddedArrayRef();
  90     ArrayRef<RVec> velocities;
  91     if (pressureCouplingType == PressureCoupling::CRescale)
  92     {
  93         velocities = statePropagatorData_->velocitiesView().unpaddedArrayRef();
  94     }
  95     // Freeze groups are not implemented
  96     ArrayRef<const unsigned short> cFreeze;
  97     // Coordinates are always scaled except for GPU update (not implemented currently)
  98     const bool scaleCoordinates = true;
  99     // Atom range
 100     const int startAtom = 0;
 101     const int numAtoms  = mdAtoms_->mdatoms()->homenr;
 102
 103     pressureCouplingScaleBoxAndCoordinates<pressureCouplingType>(
 104             inputrec_, boxScalingMatrix_, box, boxRel_, startAtom, numAtoms, positions, velocities, cFreeze, nrnb_, scaleCoordinates);
 105 }
 106
 107 void FirstOrderPressureCoupling::scheduleTask(Step step, Time /*unused*/, const RegisterRunFunction& registerRunFunction)
 108 {
 109     if (do_per_step(step + couplingFrequency_ + couplingOffset_, couplingFrequency_))
 110     {
 111         if (pressureCouplingType_ == PressureCoupling::Berendsen)
 112         {
 113             registerRunFunction([this, step]() {
 114                 calculateScalingMatrix<PressureCoupling::Berendsen>(step);
 115                 scaleBoxAndCoordinates<PressureCoupling::Berendsen>();
 116             });
 117         }
 118         else if (pressureCouplingType_ == PressureCoupling::CRescale)
 119         {
 120             registerRunFunction([this, step]() {
 121                 calculateScalingMatrix<PressureCoupling::CRescale>(step);
 122                 scaleBoxAndCoordinates<PressureCoupling::CRescale>();
 123             });
 124         }
 125     }
 126 }
 127
 128 void FirstOrderPressureCoupling::elementSetup()
 129 {
 130     if (inputrecPreserveShape(inputrec_))
 131     {
 132         auto*     box  = statePropagatorData_->box();
 133         const int ndim = inputrec_->epct == PressureCouplingType::SemiIsotropic ? 2 : 3;
 134         do_box_rel(ndim, inputrec_->deform, boxRel_, box, true);
 135     }
 136 }
 137
 138 real FirstOrderPressureCoupling::conservedEnergyContribution(Step step)
 139 {
 140     if (step == conservedEnergyContributionStep_
 141         && reportPreviousStepConservedEnergy_ == ReportPreviousStepConservedEnergy::Yes)
 142     {
 143         return previousStepConservedEnergyContribution_;
 144     }
 145     return conservedEnergyContribution_;
 146 }
 147
 148 namespace
 149 {
 150 /*!
 151  * \brief Enum describing the contents FirstOrderPressureCoupling writes to modular checkpoint
 152  *
 153  * When changing the checkpoint content, add a new element just above Count, and adjust the
 154  * checkpoint functionality.
 155  */
 156 enum class CheckpointVersion
 157 {
 158     Base, //!< First version of modular checkpointing
 159     Count //!< Number of entries. Add new versions right above this!
 160 };
 161 constexpr auto c_currentVersion = CheckpointVersion(int(CheckpointVersion::Count) - 1);
 162 } // namespace
 163
 164 template<CheckpointDataOperation operation>
 165 void FirstOrderPressureCoupling::doCheckpointData(CheckpointData<operation>* checkpointData)
 166 {
 167     checkpointVersion(checkpointData, "FirstOrderPressureCoupling version", c_currentVersion);
 168
 169     checkpointData->scalar("conserved energy contribution", &conservedEnergyContribution_);
 170     checkpointData->scalar("conserved energy step", &conservedEnergyContributionStep_);
 171     checkpointData->tensor("relative box vector", boxRel_);
 172 }
 173
 174 void FirstOrderPressureCoupling::saveCheckpointState(std::optional<WriteCheckpointData> checkpointData,
 175                                                      const t_commrec*                   cr)
 176 {
 177     if (MASTER(cr))
 178     {
 179         doCheckpointData<CheckpointDataOperation::Write>(&checkpointData.value());
 180     }
 181 }
 182
 183 void FirstOrderPressureCoupling::restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData,
 184                                                         const t_commrec*                  cr)
 185 {
 186     if (MASTER(cr))
 187     {
 188         doCheckpointData<CheckpointDataOperation::Read>(&checkpointData.value());
 189     }
 190     if (DOMAINDECOMP(cr))
 191     {
 192         dd_bcast(cr->dd, sizeof(conservedEnergyContribution_), &conservedEnergyContribution_);
 193         dd_bcast(cr->dd, sizeof(conservedEnergyContributionStep_), &conservedEnergyContributionStep_);
 194         dd_bcast(cr->dd, sizeof(boxRel_), boxRel_);
 195     }
 196 }
 197
 198 const std::string& FirstOrderPressureCoupling::clientID()
 199 {
 200     return identifier_;
 201 }
 202
 203 FirstOrderPressureCoupling::FirstOrderPressureCoupling(int                  couplingFrequency,
 204                                                        int                  couplingOffset,
 205                                                        real                 couplingTimeStep,
 206                                                        StatePropagatorData* statePropagatorData,
 207                                                        EnergyData*          energyData,
 208                                                        FILE*                fplog,
 209                                                        const t_inputrec*    inputrec,
 210                                                        const MDAtoms*       mdAtoms,
 211                                                        t_nrnb*              nrnb,
 212                                                        ReportPreviousStepConservedEnergy reportPreviousStepConservedEnergy) :
 213     pressureCouplingType_(inputrec->epc),
 214     couplingTimeStep_(couplingTimeStep),
 215     couplingFrequency_(couplingFrequency),
 216     couplingOffset_(couplingOffset),
 217     boxScalingMatrix_{ { 0 } },
 218     boxRel_{ { 0 } },
 219     conservedEnergyContribution_(0),
 220     previousStepConservedEnergyContribution_(0),
 221     conservedEnergyContributionStep_(-1),
 222     reportPreviousStepConservedEnergy_(reportPreviousStepConservedEnergy),
 223     statePropagatorData_(statePropagatorData),
 224     energyData_(energyData),
 225     fplog_(fplog),
 226     inputrec_(inputrec),
 227     mdAtoms_(mdAtoms),
 228     nrnb_(nrnb),
 229     identifier_("FirstOrderPressureCoupling-" + std::string(enumValueToString(pressureCouplingType_)))
 230 {
 231     energyData->addConservedEnergyContribution(
 232             [this](Step step, Time /*unused*/) { return conservedEnergyContribution(step); });
 233 }
 234
 235 ISimulatorElement* FirstOrderPressureCoupling::getElementPointerImpl(
 236         LegacySimulatorData*                    legacySimulatorData,
 237         ModularSimulatorAlgorithmBuilderHelper* builderHelper,
 238         StatePropagatorData*                    statePropagatorData,
 239         EnergyData*                             energyData,
 240         FreeEnergyPerturbationData gmx_unused* freeEnergyPerturbationData,
 241         GlobalCommunicationHelper gmx_unused* globalCommunicationHelper,
 242         int                                   offset,
 243         ReportPreviousStepConservedEnergy     reportPreviousStepConservedEnergy)
 244 {
 245     return builderHelper->storeElement(std::make_unique<FirstOrderPressureCoupling>(
 246             legacySimulatorData->inputrec->nstpcouple,
 247             offset,
 248             legacySimulatorData->inputrec->delta_t * legacySimulatorData->inputrec->nstpcouple,
 249             statePropagatorData,
 250             energyData,
 251             legacySimulatorData->fplog,
 252             legacySimulatorData->inputrec,
 253             legacySimulatorData->mdAtoms,
 254             legacySimulatorData->nrnb,
 255             reportPreviousStepConservedEnergy));
 256 }
 257
 258 } // namespace gmx