2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
36 * \brief Defines the Parrinello-Rahman barostat for the modular simulator
38 * \author Pascal Merz <pascal.merz@me.com>
39 * \ingroup module_modularsimulator
44 #include "parrinellorahmanbarostat.h"
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"
58 #include "energydata.h"
59 #include "modularsimulator.h"
60 #include "simulatoralgorithm.h"
61 #include "statepropagatordata.h"
66 ParrinelloRahmanBarostat::ParrinelloRahmanBarostat(int nstpcouple,
68 real couplingTimeStep,
70 StatePropagatorData* statePropagatorData,
71 EnergyData* energyData,
73 const t_inputrec* inputrec,
74 const MDAtoms* mdAtoms) :
75 nstpcouple_(nstpcouple),
77 couplingTimeStep_(couplingTimeStep),
81 boxVelocity_{ { 0 } },
82 statePropagatorData_(statePropagatorData),
83 energyData_(energyData),
88 energyData->setParrinelloRahamnBarostat(this);
91 void ParrinelloRahmanBarostat::connectWithPropagator(const PropagatorBarostatConnection& connectionData)
93 scalingTensor_ = connectionData.getViewOnPRScalingMatrix();
94 propagatorCallback_ = connectionData.getPRScalingCallback();
97 void ParrinelloRahmanBarostat::scheduleTask(Step step,
99 const RegisterRunFunction& registerRunFunction)
101 const bool scaleOnNextStep = do_per_step(step + nstpcouple_ + offset_ + 1, nstpcouple_);
102 const bool scaleOnThisStep = do_per_step(step + nstpcouple_ + offset_, nstpcouple_);
106 registerRunFunction([this]() { scaleBoxAndPositions(); });
110 registerRunFunction([this, step]() { integrateBoxVelocityEquations(step); });
111 // let propagator know that it will have to scale on next step
112 propagatorCallback_(step + 1);
116 void ParrinelloRahmanBarostat::integrateBoxVelocityEquations(Step step)
118 auto box = statePropagatorData_->constBox();
119 parrinellorahman_pcoupl(fplog_,
123 energyData_->pressure(step),
127 scalingTensor_.data(),
130 // multiply matrix by the coupling time step to avoid having the propagator needing to know about that
131 msmul(scalingTensor_.data(), couplingTimeStep_, scalingTensor_.data());
134 void ParrinelloRahmanBarostat::scaleBoxAndPositions()
136 // Propagate the box by the box velocities
137 auto box = statePropagatorData_->box();
138 for (int i = 0; i < DIM; i++)
140 for (int m = 0; m <= i; m++)
142 box[i][m] += couplingTimeStep_ * boxVelocity_[i][m];
145 preserve_box_shape(inputrec_, boxRel_, box);
147 // Scale the coordinates
149 const int homenr = mdAtoms_->mdatoms()->homenr;
150 auto x = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
151 for (int n = start; n < start + homenr; n++)
153 tmvmul_ur0(mu_, x[n], x[n]);
157 void ParrinelloRahmanBarostat::elementSetup()
159 if (!propagatorCallback_ || scalingTensor_.empty())
161 throw MissingElementConnectionError(
162 "Parrinello-Rahman barostat was not connected to a propagator.\n"
163 "Connection to a propagator element is needed to scale the velocities.\n"
164 "Use connectWithPropagator(...) before building the ModularSimulatorAlgorithm "
168 if (inputrecPreserveShape(inputrec_))
170 auto box = statePropagatorData_->box();
171 const int ndim = inputrec_->epct == PressureCouplingType::SemiIsotropic ? 2 : 3;
172 do_box_rel(ndim, inputrec_->deform, boxRel_, box, true);
175 const bool scaleOnInitStep = do_per_step(initStep_ + nstpcouple_ + offset_, nstpcouple_);
178 // If we need to scale on the first step, we need to set the scaling matrix using the current
179 // box velocity. If this is a fresh start, we will hence not move the box (this does currently
180 // never happen as the offset is set to -1 in all cases). If this is a restart, we will use
181 // the saved box velocity which we would have updated right before checkpointing.
182 // Setting bFirstStep = true in parrinellorahman_pcoupl (last argument) makes sure that only
183 // the scaling matrix is calculated, without updating the box velocities.
184 // The call to parrinellorahman_pcoupl is using nullptr for fplog (since we don't expect any
185 // output here) and for the pressure (since it might not be calculated yet, and we don't need it).
186 auto box = statePropagatorData_->constBox();
187 parrinellorahman_pcoupl(nullptr,
195 scalingTensor_.data(),
198 // multiply matrix by the coupling time step to avoid having the propagator needing to know about that
199 msmul(scalingTensor_.data(), couplingTimeStep_, scalingTensor_.data());
201 propagatorCallback_(initStep_);
205 const rvec* ParrinelloRahmanBarostat::boxVelocities() const
210 real ParrinelloRahmanBarostat::conservedEnergyContribution() const
213 const auto* box = statePropagatorData_->constBox();
214 real maxBoxLength = std::max({ box[XX][XX], box[YY][YY], box[ZZ][ZZ] });
215 real volume = det(box);
217 // contribution from the pressure momenta
218 for (int i = 0; i < DIM; i++)
220 for (int j = 0; j <= i; j++)
222 real invMass = c_presfac * (4 * M_PI * M_PI * inputrec_->compress[i][j])
223 / (3 * inputrec_->tau_p * inputrec_->tau_p * maxBoxLength);
226 energy += 0.5 * boxVelocity_[i][j] * boxVelocity_[i][j] / invMass;
231 /* Contribution from the PV term.
232 * Note that with non-zero off-diagonal reference pressures,
233 * i.e. applied shear stresses, there are additional terms.
234 * We don't support this here, since that requires keeping
235 * track of unwrapped box diagonal elements. This case is
236 * excluded in integratorHasConservedEnergyQuantity().
238 energy += volume * trace(inputrec_->ref_p) / (DIM * c_presfac);
246 * \brief Enum describing the contents ParrinelloRahmanBarostat writes to modular checkpoint
248 * When changing the checkpoint content, add a new element just above Count, and adjust the
249 * checkpoint functionality.
251 enum class CheckpointVersion
253 Base, //!< First version of modular checkpointing
254 Count //!< Number of entries. Add new versions right above this!
256 constexpr auto c_currentVersion = CheckpointVersion(int(CheckpointVersion::Count) - 1);
259 template<CheckpointDataOperation operation>
260 void ParrinelloRahmanBarostat::doCheckpointData(CheckpointData<operation>* checkpointData)
262 checkpointVersion(checkpointData, "ParrinelloRahmanBarostat version", c_currentVersion);
264 checkpointData->tensor("box velocity", boxVelocity_);
265 checkpointData->tensor("relative box vector", boxRel_);
268 void ParrinelloRahmanBarostat::saveCheckpointState(std::optional<WriteCheckpointData> checkpointData,
273 doCheckpointData<CheckpointDataOperation::Write>(&checkpointData.value());
277 void ParrinelloRahmanBarostat::restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData,
282 doCheckpointData<CheckpointDataOperation::Read>(&checkpointData.value());
284 if (DOMAINDECOMP(cr))
286 dd_bcast(cr->dd, sizeof(boxVelocity_), boxVelocity_);
287 dd_bcast(cr->dd, sizeof(boxRel_), boxRel_);
291 const std::string& ParrinelloRahmanBarostat::clientID()
296 ISimulatorElement* ParrinelloRahmanBarostat::getElementPointerImpl(
297 LegacySimulatorData* legacySimulatorData,
298 ModularSimulatorAlgorithmBuilderHelper* builderHelper,
299 StatePropagatorData* statePropagatorData,
300 EnergyData* energyData,
301 FreeEnergyPerturbationData gmx_unused* freeEnergyPerturbationData,
302 GlobalCommunicationHelper gmx_unused* globalCommunicationHelper,
305 auto* element = builderHelper->storeElement(std::make_unique<ParrinelloRahmanBarostat>(
306 legacySimulatorData->inputrec->nstpcouple,
308 legacySimulatorData->inputrec->delta_t * legacySimulatorData->inputrec->nstpcouple,
309 legacySimulatorData->inputrec->init_step,
312 legacySimulatorData->fplog,
313 legacySimulatorData->inputrec,
314 legacySimulatorData->mdAtoms));
315 auto* barostat = static_cast<ParrinelloRahmanBarostat*>(element);
316 builderHelper->registerBarostat([barostat](const PropagatorBarostatConnection& connection) {
317 barostat->connectWithPropagator(connection);