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 constraint element for the modular simulator
38 * \author Pascal Merz <pascal.merz@me.com>
39 * \ingroup module_modularsimulator
44 #include "constraintelement.h"
46 #include "gromacs/math/vec.h"
47 #include "gromacs/mdlib/mdatoms.h"
48 #include "gromacs/mdtypes/commrec.h"
49 #include "gromacs/mdtypes/enerdata.h"
50 #include "gromacs/mdtypes/inputrec.h"
51 #include "gromacs/mdtypes/mdatom.h"
52 #include "gromacs/mdtypes/state.h"
53 #include "gromacs/utility/fatalerror.h"
55 #include "energydata.h"
56 #include "freeenergyperturbationdata.h"
57 #include "modularsimulator.h"
58 #include "simulatoralgorithm.h"
59 #include "statepropagatordata.h"
63 template<ConstraintVariable variable>
64 ConstraintsElement<variable>::ConstraintsElement(Constraints* constr,
65 StatePropagatorData* statePropagatorData,
66 EnergyData* energyData,
67 FreeEnergyPerturbationData* freeEnergyPerturbationData,
70 const t_inputrec* inputrec,
71 const t_mdatoms* mdAtoms) :
72 nextVirialCalculationStep_(-1),
73 nextEnergyWritingStep_(-1),
74 nextLogWritingStep_(-1),
75 isMasterRank_(isMaster),
76 statePropagatorData_(statePropagatorData),
77 energyData_(energyData),
78 freeEnergyPerturbationData_(freeEnergyPerturbationData),
84 GMX_ASSERT(constr_, "Constraint element created but constr == nullptr");
87 template<ConstraintVariable variable>
88 void ConstraintsElement<variable>::elementSetup()
90 if (!inputrec_->bContinuation
91 && ((variable == ConstraintVariable::Positions && inputrec_->eI == IntegrationAlgorithm::MD)
92 || (variable == ConstraintVariable::Velocities && inputrec_->eI == IntegrationAlgorithm::VV)))
94 const real lambdaBonded =
95 freeEnergyPerturbationData_
96 ? freeEnergyPerturbationData_->constLambdaView()[static_cast<int>(
97 FreeEnergyPerturbationCouplingType::Bonded)]
99 // Constrain the initial coordinates and velocities
100 do_constrain_first(fplog_,
103 statePropagatorData_->totalNumAtoms(),
104 statePropagatorData_->localNumAtoms(),
105 statePropagatorData_->positionsView(),
106 statePropagatorData_->velocitiesView(),
107 statePropagatorData_->box(),
112 if (inputrec_->eConstrAlg == ConstraintAlgorithm::Lincs)
115 "RMS relative constraint deviation after constraining: %.2e\n",
122 template<ConstraintVariable variable>
123 void ConstraintsElement<variable>::scheduleTask(Step step,
124 Time gmx_unused time,
125 const RegisterRunFunction& registerRunFunction)
127 bool calculateVirial = (step == nextVirialCalculationStep_);
128 bool writeLog = (step == nextLogWritingStep_);
129 bool writeEnergy = (step == nextEnergyWritingStep_);
131 // register constraining
132 registerRunFunction([this, step, calculateVirial, writeLog, writeEnergy]() {
133 apply(step, calculateVirial, writeLog, writeEnergy);
137 template<ConstraintVariable variable>
138 void ConstraintsElement<variable>::apply(Step step, bool calculateVirial, bool writeLog, bool writeEnergy)
142 ArrayRefWithPadding<RVec> x;
143 ArrayRefWithPadding<RVec> xprime;
144 ArrayRef<RVec> min_proj;
145 ArrayRefWithPadding<RVec> v;
147 const real lambdaBonded =
148 freeEnergyPerturbationData_
149 ? freeEnergyPerturbationData_
150 ->constLambdaView()[static_cast<int>(FreeEnergyPerturbationCouplingType::Bonded)]
156 case ConstraintVariable::Positions:
157 x = statePropagatorData_->previousPositionsView();
158 xprime = statePropagatorData_->positionsView();
159 v = statePropagatorData_->velocitiesView();
161 case ConstraintVariable::Velocities:
162 x = statePropagatorData_->positionsView();
163 xprime = statePropagatorData_->velocitiesView();
164 min_proj = statePropagatorData_->velocitiesView().unpaddedArrayRef();
166 default: gmx_fatal(FARGS, "Constraint algorithm not implemented for modular simulator.");
169 constr_->apply(writeLog,
177 statePropagatorData_->box(),
187 if (inputrec_->eI == IntegrationAlgorithm::VV)
189 // For some reason, the shake virial in VV is reset twice a step.
190 // Energy element will only do this once per step.
191 // TODO: Investigate this
192 clear_mat(energyData_->constraintVirial(step));
194 energyData_->addToConstraintVirial(vir_con, step);
197 /* The factor of 2 correction is necessary because half of the constraint
198 * force is removed in the VV step. This factor is either exact or a very
199 * good approximation, statistically insignificant in any real free energy
200 * calculation. Any possible error is not a simulation propagation error,
201 * but a potential reporting error in the data that goes to dh/dlambda.
204 const real c_dvdlConstraintCorrectionFactor = EI_VV(inputrec_->eI) ? 2.0 : 1.0;
205 energyData_->enerdata()->term[F_DVDL_CONSTR] += c_dvdlConstraintCorrectionFactor * dvdlambda;
208 template<ConstraintVariable variable>
209 std::optional<SignallerCallback> ConstraintsElement<variable>::registerEnergyCallback(EnergySignallerEvent event)
211 if (event == EnergySignallerEvent::VirialCalculationStep)
213 return [this](Step step, Time /*unused*/) { nextVirialCalculationStep_ = step; };
218 template<ConstraintVariable variable>
219 std::optional<SignallerCallback> ConstraintsElement<variable>::registerTrajectorySignallerCallback(TrajectoryEvent event)
221 if (event == TrajectoryEvent::EnergyWritingStep)
223 return [this](Step step, Time /*unused*/) { nextEnergyWritingStep_ = step; };
228 template<ConstraintVariable variable>
229 std::optional<SignallerCallback> ConstraintsElement<variable>::registerLoggingCallback()
231 return [this](Step step, Time /*unused*/) { nextLogWritingStep_ = step; };
234 template<ConstraintVariable variable>
235 ISimulatorElement* ConstraintsElement<variable>::getElementPointerImpl(
236 LegacySimulatorData* legacySimulatorData,
237 ModularSimulatorAlgorithmBuilderHelper* builderHelper,
238 StatePropagatorData* statePropagatorData,
239 EnergyData* energyData,
240 FreeEnergyPerturbationData* freeEnergyPerturbationData,
241 GlobalCommunicationHelper gmx_unused* globalCommunicationHelper)
243 return builderHelper->storeElement(
244 std::make_unique<ConstraintsElement<variable>>(legacySimulatorData->constr,
247 freeEnergyPerturbationData,
248 MASTER(legacySimulatorData->cr),
249 legacySimulatorData->fplog,
250 legacySimulatorData->inputrec,
251 legacySimulatorData->mdAtoms->mdatoms()));
254 // Explicit template initializations
255 template class ConstraintsElement<ConstraintVariable::Positions>;
256 template class ConstraintsElement<ConstraintVariable::Velocities>;