src/gromacs/modularsimulator/parrinellorahmanbarostat.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020,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 Parrinello-Rahman barostat 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 "parrinellorahmanbarostat.h"
  45
  46 #include "gromacs/domdec/domdec_network.h"
  47 #include "gromacs/math/units.h"
  48 #include "gromacs/math/vec.h"
  49 #include "gromacs/mdlib/coupling.h"
  50 #include "gromacs/mdlib/mdatoms.h"
  51 #include "gromacs/mdlib/stat.h"
  52 #include "gromacs/mdtypes/checkpointdata.h"
  53 #include "gromacs/mdtypes/commrec.h"
  54 #include "gromacs/mdtypes/inputrec.h"
  55 #include "gromacs/mdtypes/mdatom.h"
  56 #include "gromacs/pbcutil/boxutilities.h"
  57
  58 #include "energydata.h"
  59 #include "modularsimulator.h"
  60 #include "simulatoralgorithm.h"
  61 #include "statepropagatordata.h"
  62
  63 namespace gmx
  64 {
  65
  66 ParrinelloRahmanBarostat::ParrinelloRahmanBarostat(int                  nstpcouple,
  67                                                    int                  offset,
  68                                                    real                 couplingTimeStep,
  69                                                    Step                 initStep,
  70                                                    StatePropagatorData* statePropagatorData,
  71                                                    EnergyData*          energyData,
  72                                                    FILE*                fplog,
  73                                                    const t_inputrec*    inputrec,
  74                                                    const MDAtoms*       mdAtoms) :
  75     nstpcouple_(nstpcouple),
  76     offset_(offset),
  77     couplingTimeStep_(couplingTimeStep),
  78     initStep_(initStep),
  79     mu_{ { 0 } },
  80     boxRel_{ { 0 } },
  81     boxVelocity_{ { 0 } },
  82     statePropagatorData_(statePropagatorData),
  83     energyData_(energyData),
  84     nextEnergyCalculationStep_(-1),
  85     fplog_(fplog),
  86     inputrec_(inputrec),
  87     mdAtoms_(mdAtoms)
  88 {
  89     energyData->setParrinelloRahmanBoxVelocities([this]() { return boxVelocity_; });
  90     energyData->addConservedEnergyContribution([this](Step gmx_used_in_debug step, Time /*unused*/) {
  91         GMX_ASSERT(conservedEnergyContributionStep_ == step,
  92                    "Parrinello-Rahman conserved energy step mismatch.");
  93         return conservedEnergyContribution_;
  94     });
  95 }
  96
  97 void ParrinelloRahmanBarostat::connectWithMatchingPropagator(const PropagatorConnection& connectionData,
  98                                                              const PropagatorTag& propagatorTag)
  99 {
 100     if (connectionData.tag == propagatorTag)
 101     {
 102         GMX_RELEASE_ASSERT(connectionData.hasParrinelloRahmanScaling(),
 103                            "Connection data lacks Parrinello-Rahman scaling");
 104         scalingTensor_      = connectionData.getViewOnPRScalingMatrix();
 105         propagatorCallback_ = connectionData.getPRScalingCallback();
 106     }
 107 }
 108
 109 void ParrinelloRahmanBarostat::scheduleTask(Step                       step,
 110                                             Time gmx_unused            time,
 111                                             const RegisterRunFunction& registerRunFunction)
 112 {
 113     const bool scaleOnNextStep = do_per_step(step + nstpcouple_ + offset_ + 1, nstpcouple_);
 114     const bool scaleOnThisStep = do_per_step(step + nstpcouple_ + offset_, nstpcouple_);
 115     const bool contributeEnergyThisStep = (step == nextEnergyCalculationStep_);
 116
 117     if (contributeEnergyThisStep)
 118     {
 119         // For compatibility with legacy md, we store this before integrating the box velocities
 120         registerRunFunction([this, step]() {
 121             conservedEnergyContribution_     = conservedEnergyContribution();
 122             conservedEnergyContributionStep_ = step;
 123         });
 124     }
 125     if (scaleOnThisStep)
 126     {
 127         registerRunFunction([this]() { scaleBoxAndPositions(); });
 128     }
 129     if (scaleOnNextStep)
 130     {
 131         registerRunFunction([this, step]() { integrateBoxVelocityEquations(step); });
 132         // let propagator know that it will have to scale on next step
 133         propagatorCallback_(step + 1);
 134     }
 135 }
 136
 137 void ParrinelloRahmanBarostat::integrateBoxVelocityEquations(Step step)
 138 {
 139     auto box = statePropagatorData_->constBox();
 140     parrinellorahman_pcoupl(fplog_,
 141                             step,
 142                             inputrec_,
 143                             couplingTimeStep_,
 144                             energyData_->pressure(step),
 145                             box,
 146                             boxRel_,
 147                             boxVelocity_,
 148                             scalingTensor_.data(),
 149                             mu_,
 150                             false);
 151     // multiply matrix by the coupling time step to avoid having the propagator needing to know about that
 152     msmul(scalingTensor_.data(), couplingTimeStep_, scalingTensor_.data());
 153 }
 154
 155 void ParrinelloRahmanBarostat::scaleBoxAndPositions()
 156 {
 157     // Propagate the box by the box velocities
 158     auto box = statePropagatorData_->box();
 159     for (int i = 0; i < DIM; i++)
 160     {
 161         for (int m = 0; m <= i; m++)
 162         {
 163             box[i][m] += couplingTimeStep_ * boxVelocity_[i][m];
 164         }
 165     }
 166     preserve_box_shape(inputrec_, boxRel_, box);
 167
 168     // Scale the coordinates
 169     const int start  = 0;
 170     const int homenr = mdAtoms_->mdatoms()->homenr;
 171     auto      x      = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
 172     for (int n = start; n < start + homenr; n++)
 173     {
 174         tmvmul_ur0(mu_, x[n], x[n]);
 175     }
 176 }
 177
 178 void ParrinelloRahmanBarostat::elementSetup()
 179 {
 180     if (!propagatorCallback_ || scalingTensor_.empty())
 181     {
 182         throw MissingElementConnectionError(
 183                 "Parrinello-Rahman barostat was not connected to a propagator.\n"
 184                 "Connection to a propagator element is needed to scale the velocities.\n"
 185                 "Use connectWithMatchingPropagator(...) before building the "
 186                 "ModularSimulatorAlgorithm "
 187                 "object.");
 188     }
 189
 190     if (inputrecPreserveShape(inputrec_))
 191     {
 192         auto      box  = statePropagatorData_->box();
 193         const int ndim = inputrec_->epct == PressureCouplingType::SemiIsotropic ? 2 : 3;
 194         do_box_rel(ndim, inputrec_->deform, boxRel_, box, true);
 195     }
 196
 197     const bool scaleOnInitStep = do_per_step(initStep_ + nstpcouple_ + offset_, nstpcouple_);
 198     if (scaleOnInitStep)
 199     {
 200         // If we need to scale on the first step, we need to set the scaling matrix using the current
 201         // box velocity. If this is a fresh start, we will hence not move the box (this does currently
 202         // never happen as the offset is set to -1 in all cases). If this is a restart, we will use
 203         // the saved box velocity which we would have updated right before checkpointing.
 204         // Setting bFirstStep = true in parrinellorahman_pcoupl (last argument) makes sure that only
 205         // the scaling matrix is calculated, without updating the box velocities.
 206         // The call to parrinellorahman_pcoupl is using nullptr for fplog (since we don't expect any
 207         // output here) and for the pressure (since it might not be calculated yet, and we don't need it).
 208         auto box = statePropagatorData_->constBox();
 209         parrinellorahman_pcoupl(nullptr,
 210                                 initStep_,
 211                                 inputrec_,
 212                                 couplingTimeStep_,
 213                                 nullptr,
 214                                 box,
 215                                 boxRel_,
 216                                 boxVelocity_,
 217                                 scalingTensor_.data(),
 218                                 mu_,
 219                                 true);
 220         // multiply matrix by the coupling time step to avoid having the propagator needing to know about that
 221         msmul(scalingTensor_.data(), couplingTimeStep_, scalingTensor_.data());
 222
 223         propagatorCallback_(initStep_);
 224     }
 225 }
 226
 227 const rvec* ParrinelloRahmanBarostat::boxVelocities() const
 228 {
 229     return boxVelocity_;
 230 }
 231
 232 real ParrinelloRahmanBarostat::conservedEnergyContribution() const
 233 {
 234     real        energy       = 0;
 235     const auto* box          = statePropagatorData_->constBox();
 236     real        maxBoxLength = std::max({ box[XX][XX], box[YY][YY], box[ZZ][ZZ] });
 237     real        volume       = det(box);
 238
 239     // contribution from the pressure momenta
 240     for (int i = 0; i < DIM; i++)
 241     {
 242         for (int j = 0; j <= i; j++)
 243         {
 244             real invMass = c_presfac * (4 * M_PI * M_PI * inputrec_->compress[i][j])
 245                            / (3 * inputrec_->tau_p * inputrec_->tau_p * maxBoxLength);
 246             if (invMass > 0)
 247             {
 248                 energy += 0.5 * boxVelocity_[i][j] * boxVelocity_[i][j] / invMass;
 249             }
 250         }
 251     }
 252
 253     /* Contribution from the PV term.
 254      * Note that with non-zero off-diagonal reference pressures,
 255      * i.e. applied shear stresses, there are additional terms.
 256      * We don't support this here, since that requires keeping
 257      * track of unwrapped box diagonal elements. This case is
 258      * excluded in integratorHasConservedEnergyQuantity().
 259      */
 260     energy += volume * trace(inputrec_->ref_p) / (DIM * c_presfac);
 261
 262     return energy;
 263 }
 264
 265 namespace
 266 {
 267 /*!
 268  * \brief Enum describing the contents ParrinelloRahmanBarostat writes to modular checkpoint
 269  *
 270  * When changing the checkpoint content, add a new element just above Count, and adjust the
 271  * checkpoint functionality.
 272  */
 273 enum class CheckpointVersion
 274 {
 275     Base, //!< First version of modular checkpointing
 276     Count //!< Number of entries. Add new versions right above this!
 277 };
 278 constexpr auto c_currentVersion = CheckpointVersion(int(CheckpointVersion::Count) - 1);
 279 } // namespace
 280
 281 template<CheckpointDataOperation operation>
 282 void ParrinelloRahmanBarostat::doCheckpointData(CheckpointData<operation>* checkpointData)
 283 {
 284     checkpointVersion(checkpointData, "ParrinelloRahmanBarostat version", c_currentVersion);
 285
 286     checkpointData->tensor("box velocity", boxVelocity_);
 287     checkpointData->tensor("relative box vector", boxRel_);
 288 }
 289
 290 void ParrinelloRahmanBarostat::saveCheckpointState(std::optional<WriteCheckpointData> checkpointData,
 291                                                    const t_commrec*                   cr)
 292 {
 293     if (MASTER(cr))
 294     {
 295         doCheckpointData<CheckpointDataOperation::Write>(&checkpointData.value());
 296     }
 297 }
 298
 299 void ParrinelloRahmanBarostat::restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData,
 300                                                       const t_commrec*                  cr)
 301 {
 302     if (MASTER(cr))
 303     {
 304         doCheckpointData<CheckpointDataOperation::Read>(&checkpointData.value());
 305     }
 306     if (DOMAINDECOMP(cr))
 307     {
 308         dd_bcast(cr->dd, sizeof(boxVelocity_), boxVelocity_);
 309         dd_bcast(cr->dd, sizeof(boxRel_), boxRel_);
 310     }
 311 }
 312
 313 const std::string& ParrinelloRahmanBarostat::clientID()
 314 {
 315     return identifier_;
 316 }
 317
 318 std::optional<SignallerCallback> ParrinelloRahmanBarostat::registerEnergyCallback(EnergySignallerEvent event)
 319 {
 320     if (event == EnergySignallerEvent::EnergyCalculationStep)
 321     {
 322         return [this](Step step, Time /*unused*/) { nextEnergyCalculationStep_ = step; };
 323     }
 324     return std::nullopt;
 325 }
 326
 327 ISimulatorElement* ParrinelloRahmanBarostat::getElementPointerImpl(
 328         LegacySimulatorData*                    legacySimulatorData,
 329         ModularSimulatorAlgorithmBuilderHelper* builderHelper,
 330         StatePropagatorData*                    statePropagatorData,
 331         EnergyData*                             energyData,
 332         FreeEnergyPerturbationData gmx_unused* freeEnergyPerturbationData,
 333         GlobalCommunicationHelper gmx_unused* globalCommunicationHelper,
 334         Offset                                offset,
 335         const PropagatorTag&                  propagatorTag)
 336 {
 337     auto* element  = builderHelper->storeElement(std::make_unique<ParrinelloRahmanBarostat>(
 338             legacySimulatorData->inputrec->nstpcouple,
 339             offset,
 340             legacySimulatorData->inputrec->delta_t * legacySimulatorData->inputrec->nstpcouple,
 341             legacySimulatorData->inputrec->init_step,
 342             statePropagatorData,
 343             energyData,
 344             legacySimulatorData->fplog,
 345             legacySimulatorData->inputrec,
 346             legacySimulatorData->mdAtoms));
 347     auto* barostat = static_cast<ParrinelloRahmanBarostat*>(element);
 348     builderHelper->registerTemperaturePressureControl(
 349             [barostat, propagatorTag](const PropagatorConnection& connection) {
 350                 barostat->connectWithMatchingPropagator(connection, propagatorTag);
 351             });
 352     return element;
 353 }
 354
 355 } // namespace gmx