2 * This file is part of the GROMACS molecular simulation package.
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.
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 element performing first-order pressure coupling for the modular simulator
38 * \author Pascal Merz <pascal.merz@me.com>
39 * \ingroup module_modularsimulator
44 #include "firstorderpressurecoupling.h"
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"
56 #include "energydata.h"
57 #include "statepropagatordata.h"
58 #include "simulatoralgorithm.h"
63 template<PressureCoupling pressureCouplingType>
64 void FirstOrderPressureCoupling::calculateScalingMatrix(Step step)
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();
71 previousStepConservedEnergyContribution_ = conservedEnergyContribution_;
72 pressureCouplingCalculateScalingMatrix<pressureCouplingType>(fplog_,
81 &conservedEnergyContribution_);
82 conservedEnergyContributionStep_ = step;
85 template<PressureCoupling pressureCouplingType>
86 void FirstOrderPressureCoupling::scaleBoxAndCoordinates()
88 auto* box = statePropagatorData_->box();
89 auto positions = statePropagatorData_->positionsView().unpaddedArrayRef();
90 ArrayRef<RVec> velocities;
91 if (pressureCouplingType == PressureCoupling::CRescale)
93 velocities = statePropagatorData_->velocitiesView().unpaddedArrayRef();
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;
100 const int startAtom = 0;
101 const int numAtoms = mdAtoms_->mdatoms()->homenr;
103 pressureCouplingScaleBoxAndCoordinates<pressureCouplingType>(
104 inputrec_, boxScalingMatrix_, box, boxRel_, startAtom, numAtoms, positions, velocities, cFreeze, nrnb_, scaleCoordinates);
107 void FirstOrderPressureCoupling::scheduleTask(Step step, Time /*unused*/, const RegisterRunFunction& registerRunFunction)
109 if (do_per_step(step + couplingFrequency_ + couplingOffset_, couplingFrequency_))
111 if (pressureCouplingType_ == PressureCoupling::Berendsen)
113 registerRunFunction([this, step]() {
114 calculateScalingMatrix<PressureCoupling::Berendsen>(step);
115 scaleBoxAndCoordinates<PressureCoupling::Berendsen>();
118 else if (pressureCouplingType_ == PressureCoupling::CRescale)
120 registerRunFunction([this, step]() {
121 calculateScalingMatrix<PressureCoupling::CRescale>(step);
122 scaleBoxAndCoordinates<PressureCoupling::CRescale>();
128 void FirstOrderPressureCoupling::elementSetup()
130 if (inputrecPreserveShape(inputrec_))
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);
138 real FirstOrderPressureCoupling::conservedEnergyContribution(Step step)
140 if (step == conservedEnergyContributionStep_
141 && reportPreviousStepConservedEnergy_ == ReportPreviousStepConservedEnergy::Yes)
143 return previousStepConservedEnergyContribution_;
145 return conservedEnergyContribution_;
151 * \brief Enum describing the contents FirstOrderPressureCoupling writes to modular checkpoint
153 * When changing the checkpoint content, add a new element just above Count, and adjust the
154 * checkpoint functionality.
156 enum class CheckpointVersion
158 Base, //!< First version of modular checkpointing
159 Count //!< Number of entries. Add new versions right above this!
161 constexpr auto c_currentVersion = CheckpointVersion(int(CheckpointVersion::Count) - 1);
164 template<CheckpointDataOperation operation>
165 void FirstOrderPressureCoupling::doCheckpointData(CheckpointData<operation>* checkpointData)
167 checkpointVersion(checkpointData, "FirstOrderPressureCoupling version", c_currentVersion);
169 checkpointData->scalar("conserved energy contribution", &conservedEnergyContribution_);
170 checkpointData->scalar("conserved energy step", &conservedEnergyContributionStep_);
171 checkpointData->tensor("relative box vector", boxRel_);
174 void FirstOrderPressureCoupling::saveCheckpointState(std::optional<WriteCheckpointData> checkpointData,
179 doCheckpointData<CheckpointDataOperation::Write>(&checkpointData.value());
183 void FirstOrderPressureCoupling::restoreCheckpointState(std::optional<ReadCheckpointData> checkpointData,
188 doCheckpointData<CheckpointDataOperation::Read>(&checkpointData.value());
190 if (haveDDAtomOrdering(*cr))
192 dd_bcast(cr->dd, sizeof(conservedEnergyContribution_), &conservedEnergyContribution_);
193 dd_bcast(cr->dd, sizeof(conservedEnergyContributionStep_), &conservedEnergyContributionStep_);
194 dd_bcast(cr->dd, sizeof(boxRel_), boxRel_);
198 const std::string& FirstOrderPressureCoupling::clientID()
203 FirstOrderPressureCoupling::FirstOrderPressureCoupling(int couplingFrequency,
205 real couplingTimeStep,
206 StatePropagatorData* statePropagatorData,
207 EnergyData* energyData,
209 const t_inputrec* inputrec,
210 const MDAtoms* mdAtoms,
212 ReportPreviousStepConservedEnergy reportPreviousStepConservedEnergy) :
213 pressureCouplingType_(inputrec->epc),
214 couplingTimeStep_(couplingTimeStep),
215 couplingFrequency_(couplingFrequency),
216 couplingOffset_(couplingOffset),
217 boxScalingMatrix_{ { 0 } },
219 conservedEnergyContribution_(0),
220 previousStepConservedEnergyContribution_(0),
221 conservedEnergyContributionStep_(-1),
222 reportPreviousStepConservedEnergy_(reportPreviousStepConservedEnergy),
223 statePropagatorData_(statePropagatorData),
224 energyData_(energyData),
229 identifier_("FirstOrderPressureCoupling-" + std::string(enumValueToString(pressureCouplingType_)))
231 energyData->addConservedEnergyContribution(
232 [this](Step step, Time /*unused*/) { return conservedEnergyContribution(step); });
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 ObservablesReducer* /*observablesReducer*/,
244 ReportPreviousStepConservedEnergy reportPreviousStepConservedEnergy)
246 return builderHelper->storeElement(std::make_unique<FirstOrderPressureCoupling>(
247 legacySimulatorData->inputrec->nstpcouple,
249 legacySimulatorData->inputrec->delta_t * legacySimulatorData->inputrec->nstpcouple,
252 legacySimulatorData->fplog,
253 legacySimulatorData->inputrec,
254 legacySimulatorData->mdAtoms,
255 legacySimulatorData->nrnb,
256 reportPreviousStepConservedEnergy));