src/gromacs/modularsimulator/parrinellorahmanbarostat.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 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     fplog_(fplog),
  85     inputrec_(inputrec),
  86     mdAtoms_(mdAtoms)
  87 {
  88     energyData->setParrinelloRahamnBarostat(this);
  89 }
  90
  91 void ParrinelloRahmanBarostat::connectWithPropagator(const PropagatorBarostatConnection& connectionData)
  92 {
  93     scalingTensor_      = connectionData.getViewOnPRScalingMatrix();
  94     propagatorCallback_ = connectionData.getPRScalingCallback();
  95 }
  96
  97 void ParrinelloRahmanBarostat::scheduleTask(Step step,
  98                                             Time gmx_unused            time,
  99                                             const RegisterRunFunction& registerRunFunction)
 100 {
 101     const bool scaleOnNextStep = do_per_step(step + nstpcouple_ + offset_ + 1, nstpcouple_);
 102     const bool scaleOnThisStep = do_per_step(step + nstpcouple_ + offset_, nstpcouple_);
 103
 104     if (scaleOnThisStep)
 105     {
 106         registerRunFunction([this]() { scaleBoxAndPositions(); });
 107     }
 108     if (scaleOnNextStep)
 109     {
 110         registerRunFunction([this, step]() { integrateBoxVelocityEquations(step); });
 111         // let propagator know that it will have to scale on next step
 112         propagatorCallback_(step + 1);
 113     }
 114 }
 115
 116 void ParrinelloRahmanBarostat::integrateBoxVelocityEquations(Step step)
 117 {
 118     auto box = statePropagatorData_->constBox();
 119     parrinellorahman_pcoupl(fplog_, step, inputrec_, couplingTimeStep_, energyData_->pressure(step),
 120                             box, boxRel_, boxVelocity_, scalingTensor_.data(), mu_, false);
 121     // multiply matrix by the coupling time step to avoid having the propagator needing to know about that
 122     msmul(scalingTensor_.data(), couplingTimeStep_, scalingTensor_.data());
 123 }
 124
 125 void ParrinelloRahmanBarostat::scaleBoxAndPositions()
 126 {
 127     // Propagate the box by the box velocities
 128     auto box = statePropagatorData_->box();
 129     for (int i = 0; i < DIM; i++)
 130     {
 131         for (int m = 0; m <= i; m++)
 132         {
 133             box[i][m] += couplingTimeStep_ * boxVelocity_[i][m];
 134         }
 135     }
 136     preserve_box_shape(inputrec_, boxRel_, box);
 137
 138     // Scale the coordinates
 139     const int start  = 0;
 140     const int homenr = mdAtoms_->mdatoms()->homenr;
 141     auto      x      = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
 142     for (int n = start; n < start + homenr; n++)
 143     {
 144         tmvmul_ur0(mu_, x[n], x[n]);
 145     }
 146 }
 147
 148 void ParrinelloRahmanBarostat::elementSetup()
 149 {
 150     if (!propagatorCallback_ || scalingTensor_.empty())
 151     {
 152         throw MissingElementConnectionError(
 153                 "Parrinello-Rahman barostat was not connected to a propagator.\n"
 154                 "Connection to a propagator element is needed to scale the velocities.\n"
 155                 "Use connectWithPropagator(...) before building the ModularSimulatorAlgorithm "
 156                 "object.");
 157     }
 158
 159     if (inputrecPreserveShape(inputrec_))
 160     {
 161         auto      box  = statePropagatorData_->box();
 162         const int ndim = inputrec_->epct == epctSEMIISOTROPIC ? 2 : 3;
 163         do_box_rel(ndim, inputrec_->deform, boxRel_, box, true);
 164     }
 165
 166     const bool scaleOnInitStep = do_per_step(initStep_ + nstpcouple_ + offset_, nstpcouple_);
 167     if (scaleOnInitStep)
 168     {
 169         // If we need to scale on the first step, we need to set the scaling matrix using the current
 170         // box velocity. If this is a fresh start, we will hence not move the box (this does currently
 171         // never happen as the offset is set to -1 in all cases). If this is a restart, we will use
 172         // the saved box velocity which we would have updated right before checkpointing.
 173         // Setting bFirstStep = true in parrinellorahman_pcoupl (last argument) makes sure that only
 174         // the scaling matrix is calculated, without updating the box velocities.
 175         // The call to parrinellorahman_pcoupl is using nullptr for fplog (since we don't expect any
 176         // output here) and for the pressure (since it might not be calculated yet, and we don't need it).
 177         auto box = statePropagatorData_->constBox();
 178         parrinellorahman_pcoupl(nullptr, initStep_, inputrec_, couplingTimeStep_, nullptr, box,
 179                                 boxRel_, boxVelocity_, scalingTensor_.data(), mu_, true);
 180         // multiply matrix by the coupling time step to avoid having the propagator needing to know about that
 181         msmul(scalingTensor_.data(), couplingTimeStep_, scalingTensor_.data());
 182
 183         propagatorCallback_(initStep_);
 184     }
 185 }
 186
 187 const rvec* ParrinelloRahmanBarostat::boxVelocities() const
 188 {
 189     return boxVelocity_;
 190 }
 191
 192 real ParrinelloRahmanBarostat::conservedEnergyContribution() const
 193 {
 194     real        energy       = 0;
 195     const auto* box          = statePropagatorData_->constBox();
 196     real        maxBoxLength = std::max({ box[XX][XX], box[YY][YY], box[ZZ][ZZ] });
 197     real        volume       = det(box);
 198
 199     // contribution from the pressure momenta
 200     for (int i = 0; i < DIM; i++)
 201     {
 202         for (int j = 0; j <= i; j++)
 203         {
 204             real invMass = PRESFAC * (4 * M_PI * M_PI * inputrec_->compress[i][j])
 205                            / (3 * inputrec_->tau_p * inputrec_->tau_p * maxBoxLength);
 206             if (invMass > 0)
 207             {
 208                 energy += 0.5 * boxVelocity_[i][j] * boxVelocity_[i][j] / invMass;
 209             }
 210         }
 211     }
 212
 213     /* Contribution from the PV term.
 214      * Note that with non-zero off-diagonal reference pressures,
 215      * i.e. applied shear stresses, there are additional terms.
 216      * We don't support this here, since that requires keeping
 217      * track of unwrapped box diagonal elements. This case is
 218      * excluded in integratorHasConservedEnergyQuantity().
 219      */
 220     energy += volume * trace(inputrec_->ref_p) / (DIM * PRESFAC);
 221
 222     return energy;
 223 }
 224
 225 namespace
 226 {
 227 /*!
 228  * \brief Enum describing the contents ParrinelloRahmanBarostat writes to modular checkpoint
 229  *
 230  * When changing the checkpoint content, add a new element just above Count, and adjust the
 231  * checkpoint functionality.
 232  */
 233 enum class CheckpointVersion
 234 {
 235     Base, //!< First version of modular checkpointing
 236     Count //!< Number of entries. Add new versions right above this!
 237 };
 238 constexpr auto c_currentVersion = CheckpointVersion(int(CheckpointVersion::Count) - 1);
 239 } // namespace
 240
 241 template<CheckpointDataOperation operation>
 242 void ParrinelloRahmanBarostat::doCheckpointData(CheckpointData<operation>* checkpointData)
 243 {
 244     checkpointVersion(checkpointData, "ParrinelloRahmanBarostat version", c_currentVersion);
 245
 246     checkpointData->tensor("box velocity", boxVelocity_);
 247     checkpointData->tensor("relative box vector", boxRel_);
 248 }
 249
 250 void ParrinelloRahmanBarostat::saveCheckpointState(std::optional<WriteCheckpointData> checkpointData,
 251                                                    const t_commrec*                   cr)
 252 {
 253     if (MASTER(cr))
 254     {
 255         doCheckpointData<CheckpointDataOperation::Write>(&checkpointData.value());
 256     }
 257 }
 258
 259 void ParrinelloRahmanBarostat::restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData,
 260                                                       const t_commrec*                  cr)
 261 {
 262     if (MASTER(cr))
 263     {
 264         doCheckpointData<CheckpointDataOperation::Read>(&checkpointData.value());
 265     }
 266     if (DOMAINDECOMP(cr))
 267     {
 268         dd_bcast(cr->dd, sizeof(boxVelocity_), boxVelocity_);
 269         dd_bcast(cr->dd, sizeof(boxRel_), boxRel_);
 270     }
 271 }
 272
 273 const std::string& ParrinelloRahmanBarostat::clientID()
 274 {
 275     return identifier_;
 276 }
 277
 278 ISimulatorElement* ParrinelloRahmanBarostat::getElementPointerImpl(
 279         LegacySimulatorData*                    legacySimulatorData,
 280         ModularSimulatorAlgorithmBuilderHelper* builderHelper,
 281         StatePropagatorData*                    statePropagatorData,
 282         EnergyData*                             energyData,
 283         FreeEnergyPerturbationData gmx_unused* freeEnergyPerturbationData,
 284         GlobalCommunicationHelper gmx_unused* globalCommunicationHelper,
 285         int                                   offset)
 286 {
 287     auto* element  = builderHelper->storeElement(std::make_unique<ParrinelloRahmanBarostat>(
 288             legacySimulatorData->inputrec->nstpcouple, offset,
 289             legacySimulatorData->inputrec->delta_t * legacySimulatorData->inputrec->nstpcouple,
 290             legacySimulatorData->inputrec->init_step, statePropagatorData, energyData,
 291             legacySimulatorData->fplog, legacySimulatorData->inputrec, legacySimulatorData->mdAtoms));
 292     auto* barostat = static_cast<ParrinelloRahmanBarostat*>(element);
 293     builderHelper->registerBarostat([barostat](const PropagatorBarostatConnection& connection) {
 294         barostat->connectWithPropagator(connection);
 295     });
 296     return element;
 297 }
 298
 299 } // namespace gmx