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 propagator element for the modular simulator
38 * \author Pascal Merz <pascal.merz@me.com>
39 * \ingroup module_modularsimulator
44 #include "propagator.h"
46 #include "gromacs/utility.h"
47 #include "gromacs/math/vec.h"
48 #include "gromacs/math/vectypes.h"
49 #include "gromacs/mdlib/gmx_omp_nthreads.h"
50 #include "gromacs/mdlib/mdatoms.h"
51 #include "gromacs/mdlib/update.h"
52 #include "gromacs/mdtypes/inputrec.h"
53 #include "gromacs/mdtypes/mdatom.h"
54 #include "gromacs/timing/wallcycle.h"
56 #include "modularsimulator.h"
57 #include "simulatoralgorithm.h"
58 #include "statepropagatordata.h"
64 // Names of integration steps, only used locally for error messages
65 constexpr EnumerationArray<IntegrationStage, const char*> integrationStepNames = {
66 "IntegrationStage::PositionsOnly", "IntegrationStage::VelocitiesOnly",
67 "IntegrationStage::LeapFrog", "IntegrationStage::VelocityVerletPositionsAndVelocities",
68 "IntegrationStage::ScaleVelocities", "IntegrationStage::ScalePositions"
73 template<NumVelocityScalingValues numStartVelocityScalingValues,
74 ParrinelloRahmanVelocityScaling parrinelloRahmanVelocityScaling,
75 NumVelocityScalingValues numEndVelocityScalingValues>
76 static void inline updateVelocities(int a,
80 const rvec* gmx_restrict invMassPerDim,
82 const rvec* gmx_restrict f,
84 const matrix matrixPR)
86 for (int d = 0; d < DIM; d++)
88 // TODO: Extract this into policy classes
89 if (numStartVelocityScalingValues != NumVelocityScalingValues::None
90 && parrinelloRahmanVelocityScaling == ParrinelloRahmanVelocityScaling::No)
92 v[a][d] *= lambdaStart;
94 if (numStartVelocityScalingValues != NumVelocityScalingValues::None
95 && parrinelloRahmanVelocityScaling == ParrinelloRahmanVelocityScaling::Diagonal)
97 v[a][d] *= (lambdaStart - diagPR[d]);
99 if (numStartVelocityScalingValues != NumVelocityScalingValues::None
100 && parrinelloRahmanVelocityScaling == ParrinelloRahmanVelocityScaling::Full)
102 v[a][d] = lambdaStart * v[a][d] - iprod(matrixPR[d], v[a]);
104 if (numStartVelocityScalingValues == NumVelocityScalingValues::None
105 && parrinelloRahmanVelocityScaling == ParrinelloRahmanVelocityScaling::Diagonal)
107 v[a][d] *= (1 - diagPR[d]);
109 if (numStartVelocityScalingValues == NumVelocityScalingValues::None
110 && parrinelloRahmanVelocityScaling == ParrinelloRahmanVelocityScaling::Full)
112 v[a][d] -= iprod(matrixPR[d], v[a]);
114 v[a][d] += f[a][d] * invMassPerDim[a][d] * dt;
115 if (numEndVelocityScalingValues != NumVelocityScalingValues::None)
117 v[a][d] *= lambdaEnd;
123 static void inline updatePositions(int a,
125 const rvec* gmx_restrict x,
126 rvec* gmx_restrict xprime,
127 const rvec* gmx_restrict v)
129 for (int d = 0; d < DIM; d++)
131 xprime[a][d] = x[a][d] + v[a][d] * dt;
136 template<NumVelocityScalingValues numStartVelocityScalingValues>
137 static void inline scaleVelocities(int a, real lambda, rvec* gmx_restrict v)
139 if (numStartVelocityScalingValues != NumVelocityScalingValues::None)
141 for (int d = 0; d < DIM; d++)
149 template<NumPositionScalingValues numPositionScalingValues>
150 static void inline scalePositions(int a, real lambda, rvec* gmx_restrict x)
152 if (numPositionScalingValues != NumPositionScalingValues::None)
154 for (int d = 0; d < DIM; d++)
161 //! Helper function diagonalizing the PR matrix if possible
162 template<ParrinelloRahmanVelocityScaling parrinelloRahmanVelocityScaling>
163 static inline bool diagonalizePRMatrix(matrix matrixPR, rvec diagPR)
165 if (parrinelloRahmanVelocityScaling != ParrinelloRahmanVelocityScaling::Full)
171 if (matrixPR[YY][XX] == 0 && matrixPR[ZZ][XX] == 0 && matrixPR[ZZ][YY] == 0)
173 diagPR[XX] = matrixPR[XX][XX];
174 diagPR[YY] = matrixPR[YY][YY];
175 diagPR[ZZ] = matrixPR[ZZ][ZZ];
185 //! Propagation (position only)
187 template<NumVelocityScalingValues numStartVelocityScalingValues,
188 ParrinelloRahmanVelocityScaling parrinelloRahmanVelocityScaling,
189 NumVelocityScalingValues numEndVelocityScalingValues,
190 NumPositionScalingValues numPositionScalingValues>
191 void Propagator<IntegrationStage::PositionsOnly>::run()
193 wallcycle_start(wcycle_, WallCycleCounter::Update);
195 auto* xp = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
196 const auto* x = as_rvec_array(statePropagatorData_->constPositionsView().paddedArrayRef().data());
197 const auto* v = as_rvec_array(statePropagatorData_->constVelocitiesView().paddedArrayRef().data());
199 int nth = gmx_omp_nthreads_get(ModuleMultiThread::Update);
200 int homenr = mdAtoms_->mdatoms()->homenr;
202 #pragma omp parallel for num_threads(nth) schedule(static) default(none) shared(nth, homenr, x, xp, v)
203 for (int th = 0; th < nth; th++)
207 int start_th, end_th;
208 getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
210 for (int a = start_th; a < end_th; a++)
212 updatePositions(a, timestep_, x, xp, v);
215 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
217 wallcycle_stop(wcycle_, WallCycleCounter::Update);
220 //! Propagation (scale position only)
222 template<NumVelocityScalingValues numStartVelocityScalingValues,
223 ParrinelloRahmanVelocityScaling parrinelloRahmanVelocityScaling,
224 NumVelocityScalingValues numEndVelocityScalingValues,
225 NumPositionScalingValues numPositionScalingValues>
226 void Propagator<IntegrationStage::ScalePositions>::run()
228 wallcycle_start(wcycle_, WallCycleCounter::Update);
230 auto* x = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
233 (numPositionScalingValues == NumPositionScalingValues::Single) ? positionScaling_[0] : 1.0;
235 int nth = gmx_omp_nthreads_get(ModuleMultiThread::Update);
236 int homenr = mdAtoms_->mdatoms()->homenr;
238 #pragma omp parallel for num_threads(nth) schedule(static) default(none) shared(nth, homenr, x) \
240 for (int th = 0; th < nth; th++)
244 int start_th, end_th;
245 getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
247 for (int a = start_th; a < end_th; a++)
249 scalePositions<numPositionScalingValues>(
251 (numPositionScalingValues == NumPositionScalingValues::Multiple)
252 ? positionScaling_[mdAtoms_->mdatoms()->cTC[a]]
257 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
259 wallcycle_stop(wcycle_, WallCycleCounter::Update);
262 //! Propagation (velocity only)
264 template<NumVelocityScalingValues numStartVelocityScalingValues,
265 ParrinelloRahmanVelocityScaling parrinelloRahmanVelocityScaling,
266 NumVelocityScalingValues numEndVelocityScalingValues,
267 NumPositionScalingValues numPositionScalingValues>
268 void Propagator<IntegrationStage::VelocitiesOnly>::run()
270 wallcycle_start(wcycle_, WallCycleCounter::Update);
272 auto* v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
273 const auto* f = as_rvec_array(statePropagatorData_->constForcesView().force().data());
274 const auto* invMassPerDim = mdAtoms_->mdatoms()->invMassPerDim;
276 const real lambdaStart = (numStartVelocityScalingValues == NumVelocityScalingValues::Single)
277 ? startVelocityScaling_[0]
279 const real lambdaEnd = (numEndVelocityScalingValues == NumVelocityScalingValues::Single)
280 ? endVelocityScaling_[0]
283 const bool isFullScalingMatrixDiagonal =
284 diagonalizePRMatrix<parrinelloRahmanVelocityScaling>(matrixPR_, diagPR_);
286 const int nth = gmx_omp_nthreads_get(ModuleMultiThread::Update);
287 const int homenr = mdAtoms_->mdatoms()->homenr;
289 // const variables could be shared, but gcc-8 & gcc-9 don't agree how to write that...
290 // https://www.gnu.org/software/gcc/gcc-9/porting_to.html -> OpenMP data sharing
291 #pragma omp parallel for num_threads(nth) schedule(static) default(none) shared(v, f, invMassPerDim) \
292 firstprivate(nth, homenr, lambdaStart, lambdaEnd, isFullScalingMatrixDiagonal)
293 for (int th = 0; th < nth; th++)
297 int start_th, end_th;
298 getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
300 for (int a = start_th; a < end_th; a++)
302 if (isFullScalingMatrixDiagonal)
304 updateVelocities<numStartVelocityScalingValues, ParrinelloRahmanVelocityScaling::Diagonal, numEndVelocityScalingValues>(
307 numStartVelocityScalingValues == NumVelocityScalingValues::Multiple
308 ? startVelocityScaling_[mdAtoms_->mdatoms()->cTC[a]]
310 numEndVelocityScalingValues == NumVelocityScalingValues::Multiple
311 ? endVelocityScaling_[mdAtoms_->mdatoms()->cTC[a]]
321 updateVelocities<numStartVelocityScalingValues, parrinelloRahmanVelocityScaling, numEndVelocityScalingValues>(
324 numStartVelocityScalingValues == NumVelocityScalingValues::Multiple
325 ? startVelocityScaling_[mdAtoms_->mdatoms()->cTC[a]]
327 numEndVelocityScalingValues == NumVelocityScalingValues::Multiple
328 ? endVelocityScaling_[mdAtoms_->mdatoms()->cTC[a]]
338 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
340 wallcycle_stop(wcycle_, WallCycleCounter::Update);
343 //! Propagation (leapfrog case - position and velocity)
345 template<NumVelocityScalingValues numStartVelocityScalingValues,
346 ParrinelloRahmanVelocityScaling parrinelloRahmanVelocityScaling,
347 NumVelocityScalingValues numEndVelocityScalingValues,
348 NumPositionScalingValues numPositionScalingValues>
349 void Propagator<IntegrationStage::LeapFrog>::run()
351 wallcycle_start(wcycle_, WallCycleCounter::Update);
353 auto* xp = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
354 const auto* x = as_rvec_array(statePropagatorData_->constPositionsView().paddedArrayRef().data());
355 auto* v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
356 const auto* f = as_rvec_array(statePropagatorData_->constForcesView().force().data());
357 const auto* invMassPerDim = mdAtoms_->mdatoms()->invMassPerDim;
359 const real lambdaStart = (numStartVelocityScalingValues == NumVelocityScalingValues::Single)
360 ? startVelocityScaling_[0]
362 const real lambdaEnd = (numEndVelocityScalingValues == NumVelocityScalingValues::Single)
363 ? endVelocityScaling_[0]
366 const bool isFullScalingMatrixDiagonal =
367 diagonalizePRMatrix<parrinelloRahmanVelocityScaling>(matrixPR_, diagPR_);
369 const int nth = gmx_omp_nthreads_get(ModuleMultiThread::Update);
370 const int homenr = mdAtoms_->mdatoms()->homenr;
372 // const variables could be shared, but gcc-8 & gcc-9 don't agree how to write that...
373 // https://www.gnu.org/software/gcc/gcc-9/porting_to.html -> OpenMP data sharing
374 #pragma omp parallel for num_threads(nth) schedule(static) default(none) \
375 shared(x, xp, v, f, invMassPerDim) \
376 firstprivate(nth, homenr, lambdaStart, lambdaEnd, isFullScalingMatrixDiagonal)
377 for (int th = 0; th < nth; th++)
381 int start_th, end_th;
382 getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
384 for (int a = start_th; a < end_th; a++)
386 if (isFullScalingMatrixDiagonal)
388 updateVelocities<numStartVelocityScalingValues, ParrinelloRahmanVelocityScaling::Diagonal, numEndVelocityScalingValues>(
391 numStartVelocityScalingValues == NumVelocityScalingValues::Multiple
392 ? startVelocityScaling_[mdAtoms_->mdatoms()->cTC[a]]
394 numEndVelocityScalingValues == NumVelocityScalingValues::Multiple
395 ? endVelocityScaling_[mdAtoms_->mdatoms()->cTC[a]]
405 updateVelocities<numStartVelocityScalingValues, parrinelloRahmanVelocityScaling, numEndVelocityScalingValues>(
408 numStartVelocityScalingValues == NumVelocityScalingValues::Multiple
409 ? startVelocityScaling_[mdAtoms_->mdatoms()->cTC[a]]
411 numEndVelocityScalingValues == NumVelocityScalingValues::Multiple
412 ? endVelocityScaling_[mdAtoms_->mdatoms()->cTC[a]]
420 updatePositions(a, timestep_, x, xp, v);
423 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
425 wallcycle_stop(wcycle_, WallCycleCounter::Update);
428 //! Propagation (velocity verlet stage 2 - velocity and position)
430 template<NumVelocityScalingValues numStartVelocityScalingValues,
431 ParrinelloRahmanVelocityScaling parrinelloRahmanVelocityScaling,
432 NumVelocityScalingValues numEndVelocityScalingValues,
433 NumPositionScalingValues numPositionScalingValues>
434 void Propagator<IntegrationStage::VelocityVerletPositionsAndVelocities>::run()
436 wallcycle_start(wcycle_, WallCycleCounter::Update);
438 auto* xp = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
439 const auto* x = as_rvec_array(statePropagatorData_->constPositionsView().paddedArrayRef().data());
440 auto* v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
441 const auto* f = as_rvec_array(statePropagatorData_->constForcesView().force().data());
442 const auto* invMassPerDim = mdAtoms_->mdatoms()->invMassPerDim;
444 const real lambdaStart = (numStartVelocityScalingValues == NumVelocityScalingValues::Single)
445 ? startVelocityScaling_[0]
447 const real lambdaEnd = (numEndVelocityScalingValues == NumVelocityScalingValues::Single)
448 ? endVelocityScaling_[0]
451 const bool isFullScalingMatrixDiagonal =
452 diagonalizePRMatrix<parrinelloRahmanVelocityScaling>(matrixPR_, diagPR_);
454 const int nth = gmx_omp_nthreads_get(ModuleMultiThread::Update);
455 const int homenr = mdAtoms_->mdatoms()->homenr;
457 // const variables could be shared, but gcc-8 & gcc-9 don't agree how to write that...
458 // https://www.gnu.org/software/gcc/gcc-9/porting_to.html -> OpenMP data sharing
459 #pragma omp parallel for num_threads(nth) schedule(static) default(none) \
460 shared(x, xp, v, f, invMassPerDim) \
461 firstprivate(nth, homenr, lambdaStart, lambdaEnd, isFullScalingMatrixDiagonal)
462 for (int th = 0; th < nth; th++)
466 int start_th, end_th;
467 getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
469 for (int a = start_th; a < end_th; a++)
471 if (isFullScalingMatrixDiagonal)
473 updateVelocities<numStartVelocityScalingValues, ParrinelloRahmanVelocityScaling::Diagonal, numEndVelocityScalingValues>(
476 numStartVelocityScalingValues == NumVelocityScalingValues::Multiple
477 ? startVelocityScaling_[mdAtoms_->mdatoms()->cTC[a]]
479 numEndVelocityScalingValues == NumVelocityScalingValues::Multiple
480 ? endVelocityScaling_[mdAtoms_->mdatoms()->cTC[a]]
490 updateVelocities<numStartVelocityScalingValues, parrinelloRahmanVelocityScaling, numEndVelocityScalingValues>(
493 numStartVelocityScalingValues == NumVelocityScalingValues::Multiple
494 ? startVelocityScaling_[mdAtoms_->mdatoms()->cTC[a]]
496 numEndVelocityScalingValues == NumVelocityScalingValues::Multiple
497 ? endVelocityScaling_[mdAtoms_->mdatoms()->cTC[a]]
505 updatePositions(a, timestep_, x, xp, v);
508 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
510 wallcycle_stop(wcycle_, WallCycleCounter::Update);
513 //! Scaling (velocity scaling only)
515 template<NumVelocityScalingValues numStartVelocityScalingValues,
516 ParrinelloRahmanVelocityScaling parrinelloRahmanVelocityScaling,
517 NumVelocityScalingValues numEndVelocityScalingValues,
518 NumPositionScalingValues numPositionScalingValues>
519 void Propagator<IntegrationStage::ScaleVelocities>::run()
521 if (numStartVelocityScalingValues == NumVelocityScalingValues::None)
525 wallcycle_start(wcycle_, WallCycleCounter::Update);
527 auto* v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
529 const real lambdaStart = (numStartVelocityScalingValues == NumVelocityScalingValues::Single)
530 ? startVelocityScaling_[0]
533 const int nth = gmx_omp_nthreads_get(ModuleMultiThread::Update);
534 const int homenr = mdAtoms_->mdatoms()->homenr;
536 // const variables could be shared, but gcc-8 & gcc-9 don't agree how to write that...
537 // https://www.gnu.org/software/gcc/gcc-9/porting_to.html -> OpenMP data sharing
538 #pragma omp parallel for num_threads(nth) schedule(static) default(none) shared(v) \
539 firstprivate(nth, homenr, lambdaStart)
540 for (int th = 0; th < nth; th++)
546 getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
548 for (int a = start_th; a < end_th; a++)
550 scaleVelocities<numStartVelocityScalingValues>(
552 numStartVelocityScalingValues == NumVelocityScalingValues::Multiple
553 ? startVelocityScaling_[mdAtoms_->mdatoms()->cTC[a]]
558 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
560 wallcycle_stop(wcycle_, WallCycleCounter::Update);
563 template<IntegrationStage integrationStage>
564 Propagator<integrationStage>::Propagator(double timestep,
565 StatePropagatorData* statePropagatorData,
566 const MDAtoms* mdAtoms,
567 gmx_wallcycle* wcycle) :
569 statePropagatorData_(statePropagatorData),
570 doSingleStartVelocityScaling_(false),
571 doGroupStartVelocityScaling_(false),
572 doSingleEndVelocityScaling_(false),
573 doGroupEndVelocityScaling_(false),
574 scalingStepVelocity_(-1),
583 template<IntegrationStage integrationStage>
584 void Propagator<integrationStage>::scheduleTask(Step step,
585 Time gmx_unused time,
586 const RegisterRunFunction& registerRunFunction)
588 const bool doSingleVScalingThisStep =
589 (doSingleStartVelocityScaling_ && (step == scalingStepVelocity_));
590 const bool doGroupVScalingThisStep = (doGroupStartVelocityScaling_ && (step == scalingStepVelocity_));
592 if (integrationStage == IntegrationStage::ScaleVelocities)
594 // IntegrationStage::ScaleVelocities only needs to run if some kind of
595 // velocity scaling is needed on the current step.
596 if (!doSingleVScalingThisStep && !doGroupVScalingThisStep)
602 if (integrationStage == IntegrationStage::ScalePositions)
604 // IntegrationStage::ScalePositions only needs to run if
605 // position scaling is needed on the current step.
606 if (step != scalingStepPosition_)
610 // Since IntegrationStage::ScalePositions is the only stage for which position scaling
611 // is implemented we handle it here to avoid enlarging the decision tree below.
612 if (doSinglePositionScaling_)
614 registerRunFunction([this]() {
615 run<NumVelocityScalingValues::None,
616 ParrinelloRahmanVelocityScaling::No,
617 NumVelocityScalingValues::None,
618 NumPositionScalingValues::Single>();
621 else if (doGroupPositionScaling_)
623 registerRunFunction([this]() {
624 run<NumVelocityScalingValues::None,
625 ParrinelloRahmanVelocityScaling::No,
626 NumVelocityScalingValues::None,
627 NumPositionScalingValues::Multiple>();
632 const bool doParrinelloRahmanThisStep = (step == scalingStepPR_);
634 if (doSingleVScalingThisStep)
636 if (doParrinelloRahmanThisStep)
638 if (doSingleEndVelocityScaling_)
640 registerRunFunction([this]() {
641 run<NumVelocityScalingValues::Single,
642 ParrinelloRahmanVelocityScaling::Full,
643 NumVelocityScalingValues::Single,
644 NumPositionScalingValues::None>();
649 registerRunFunction([this]() {
650 run<NumVelocityScalingValues::Single,
651 ParrinelloRahmanVelocityScaling::Full,
652 NumVelocityScalingValues::None,
653 NumPositionScalingValues::None>();
659 if (doSingleEndVelocityScaling_)
661 registerRunFunction([this]() {
662 run<NumVelocityScalingValues::Single,
663 ParrinelloRahmanVelocityScaling::No,
664 NumVelocityScalingValues::Single,
665 NumPositionScalingValues::None>();
670 registerRunFunction([this]() {
671 run<NumVelocityScalingValues::Single,
672 ParrinelloRahmanVelocityScaling::No,
673 NumVelocityScalingValues::None,
674 NumPositionScalingValues::None>();
679 else if (doGroupVScalingThisStep)
681 if (doParrinelloRahmanThisStep)
683 if (doGroupEndVelocityScaling_)
685 registerRunFunction([this]() {
686 run<NumVelocityScalingValues::Multiple,
687 ParrinelloRahmanVelocityScaling::Full,
688 NumVelocityScalingValues::Multiple,
689 NumPositionScalingValues::None>();
694 registerRunFunction([this]() {
695 run<NumVelocityScalingValues::Multiple,
696 ParrinelloRahmanVelocityScaling::Full,
697 NumVelocityScalingValues::None,
698 NumPositionScalingValues::None>();
704 if (doGroupEndVelocityScaling_)
706 registerRunFunction([this]() {
707 run<NumVelocityScalingValues::Multiple,
708 ParrinelloRahmanVelocityScaling::No,
709 NumVelocityScalingValues::Multiple,
710 NumPositionScalingValues::None>();
715 registerRunFunction([this]() {
716 run<NumVelocityScalingValues::Multiple,
717 ParrinelloRahmanVelocityScaling::No,
718 NumVelocityScalingValues::None,
719 NumPositionScalingValues::None>();
726 if (doParrinelloRahmanThisStep)
728 registerRunFunction([this]() {
729 run<NumVelocityScalingValues::None,
730 ParrinelloRahmanVelocityScaling::Full,
731 NumVelocityScalingValues::None,
732 NumPositionScalingValues::None>();
737 registerRunFunction([this]() {
738 run<NumVelocityScalingValues::None,
739 ParrinelloRahmanVelocityScaling::No,
740 NumVelocityScalingValues::None,
741 NumPositionScalingValues::None>();
747 template<IntegrationStage integrationStage>
748 constexpr bool hasStartVelocityScaling()
750 return (integrationStage == IntegrationStage::VelocitiesOnly
751 || integrationStage == IntegrationStage::LeapFrog
752 || integrationStage == IntegrationStage::VelocityVerletPositionsAndVelocities
753 || integrationStage == IntegrationStage::ScaleVelocities);
756 template<IntegrationStage integrationStage>
757 constexpr bool hasEndVelocityScaling()
759 return (hasStartVelocityScaling<integrationStage>()
760 && integrationStage != IntegrationStage::ScaleVelocities);
763 template<IntegrationStage integrationStage>
764 constexpr bool hasPositionScaling()
766 return (integrationStage == IntegrationStage::ScalePositions);
769 template<IntegrationStage integrationStage>
770 constexpr bool hasParrinelloRahmanScaling()
772 return (integrationStage == IntegrationStage::VelocitiesOnly
773 || integrationStage == IntegrationStage::LeapFrog
774 || integrationStage == IntegrationStage::VelocityVerletPositionsAndVelocities);
777 template<IntegrationStage integrationStage>
778 void Propagator<integrationStage>::setNumVelocityScalingVariables(int numVelocityScalingVariables,
779 ScaleVelocities scaleVelocities)
782 hasStartVelocityScaling<integrationStage>() || hasEndVelocityScaling<integrationStage>(),
783 formatString("Velocity scaling not implemented for %s", integrationStepNames[integrationStage])
785 GMX_RELEASE_ASSERT(startVelocityScaling_.empty(),
786 "Number of velocity scaling variables cannot be changed once set.");
788 const bool scaleEndVelocities = (scaleVelocities == ScaleVelocities::PreStepAndPostStep);
789 startVelocityScaling_.resize(numVelocityScalingVariables, 1.);
790 if (scaleEndVelocities)
792 endVelocityScaling_.resize(numVelocityScalingVariables, 1.);
794 doSingleStartVelocityScaling_ = numVelocityScalingVariables == 1;
795 doGroupStartVelocityScaling_ = numVelocityScalingVariables > 1;
796 doSingleEndVelocityScaling_ = doSingleStartVelocityScaling_ && scaleEndVelocities;
797 doGroupEndVelocityScaling_ = doGroupStartVelocityScaling_ && scaleEndVelocities;
800 template<IntegrationStage integrationStage>
801 void Propagator<integrationStage>::setNumPositionScalingVariables(int numPositionScalingVariables)
803 GMX_RELEASE_ASSERT(hasPositionScaling<integrationStage>(),
804 formatString("Position scaling not implemented for %s",
805 integrationStepNames[integrationStage])
807 GMX_RELEASE_ASSERT(positionScaling_.empty(),
808 "Number of position scaling variables cannot be changed once set.");
809 positionScaling_.resize(numPositionScalingVariables, 1.);
810 doSinglePositionScaling_ = (numPositionScalingVariables == 1);
811 doGroupPositionScaling_ = (numPositionScalingVariables > 1);
814 template<IntegrationStage integrationStage>
815 ArrayRef<real> Propagator<integrationStage>::viewOnStartVelocityScaling()
817 GMX_RELEASE_ASSERT(hasStartVelocityScaling<integrationStage>(),
818 formatString("Start velocity scaling not implemented for %s",
819 integrationStepNames[integrationStage])
821 GMX_RELEASE_ASSERT(!startVelocityScaling_.empty(),
822 "Number of velocity scaling variables not set.");
824 return startVelocityScaling_;
827 template<IntegrationStage integrationStage>
828 ArrayRef<real> Propagator<integrationStage>::viewOnEndVelocityScaling()
830 GMX_RELEASE_ASSERT(hasEndVelocityScaling<integrationStage>(),
831 formatString("End velocity scaling not implemented for %s",
832 integrationStepNames[integrationStage])
834 GMX_RELEASE_ASSERT(!endVelocityScaling_.empty(),
835 "Number of velocity scaling variables not set.");
837 return endVelocityScaling_;
840 template<IntegrationStage integrationStage>
841 ArrayRef<real> Propagator<integrationStage>::viewOnPositionScaling()
843 GMX_RELEASE_ASSERT(hasPositionScaling<integrationStage>(),
844 formatString("Position scaling not implemented for %s",
845 integrationStepNames[integrationStage])
847 GMX_RELEASE_ASSERT(!positionScaling_.empty(), "Number of position scaling variables not set.");
849 return positionScaling_;
852 template<IntegrationStage integrationStage>
853 PropagatorCallback Propagator<integrationStage>::velocityScalingCallback()
856 hasStartVelocityScaling<integrationStage>() || hasEndVelocityScaling<integrationStage>(),
857 formatString("Velocity scaling not implemented for %s", integrationStepNames[integrationStage])
860 return [this](Step step) { scalingStepVelocity_ = step; };
863 template<IntegrationStage integrationStage>
864 PropagatorCallback Propagator<integrationStage>::positionScalingCallback()
866 GMX_RELEASE_ASSERT(hasPositionScaling<integrationStage>(),
867 formatString("Position scaling not implemented for %s",
868 integrationStepNames[integrationStage])
871 return [this](Step step) { scalingStepPosition_ = step; };
874 template<IntegrationStage integrationStage>
875 ArrayRef<rvec> Propagator<integrationStage>::viewOnPRScalingMatrix()
877 GMX_RELEASE_ASSERT(hasParrinelloRahmanScaling<integrationStage>(),
878 formatString("Parrinello-Rahman scaling not implemented for %s",
879 integrationStepNames[integrationStage])
882 clear_mat(matrixPR_);
883 // gcc-5 needs this to be explicit (all other tested compilers would be ok
884 // with simply returning matrixPR)
885 return ArrayRef<rvec>(matrixPR_);
888 template<IntegrationStage integrationStage>
889 PropagatorCallback Propagator<integrationStage>::prScalingCallback()
891 GMX_RELEASE_ASSERT(hasParrinelloRahmanScaling<integrationStage>(),
892 formatString("Parrinello-Rahman scaling not implemented for %s",
893 integrationStepNames[integrationStage])
896 return [this](Step step) { scalingStepPR_ = step; };
899 template<IntegrationStage integrationStage>
900 static PropagatorConnection getConnection(Propagator<integrationStage> gmx_unused* propagator,
901 const PropagatorTag& propagatorTag)
903 // gmx_unused is needed because gcc-7 & gcc-8 can't see that
904 // propagator is used for all IntegrationStage options
906 PropagatorConnection propagatorConnection{ propagatorTag };
908 // The clang-tidy version on our current CI throws 3 different warnings
909 // for the if constexpr lines, so disable linting for now. Also, this only
910 // works if the brace is on the same line, so turn off clang-format as well
913 if constexpr (hasStartVelocityScaling<integrationStage>() || hasEndVelocityScaling<integrationStage>()) {
915 propagatorConnection.setNumVelocityScalingVariables =
916 [propagator](int num, ScaleVelocities scaleVelocities) {
917 propagator->setNumVelocityScalingVariables(num, scaleVelocities);
919 propagatorConnection.getVelocityScalingCallback = [propagator]() {
920 return propagator->velocityScalingCallback();
925 if constexpr (hasStartVelocityScaling<integrationStage>()) {
927 propagatorConnection.getViewOnStartVelocityScaling = [propagator]() {
928 return propagator->viewOnStartVelocityScaling();
933 if constexpr (hasEndVelocityScaling<integrationStage>()) {
935 propagatorConnection.getViewOnEndVelocityScaling = [propagator]() {
936 return propagator->viewOnEndVelocityScaling();
941 if constexpr (hasPositionScaling<integrationStage>()) {
943 propagatorConnection.setNumPositionScalingVariables = [propagator](int num) {
944 propagator->setNumPositionScalingVariables(num);
946 propagatorConnection.getViewOnPositionScaling = [propagator]() {
947 return propagator->viewOnPositionScaling();
949 propagatorConnection.getPositionScalingCallback = [propagator]() {
950 return propagator->positionScalingCallback();
955 if constexpr (hasParrinelloRahmanScaling<integrationStage>()) {
957 propagatorConnection.getViewOnPRScalingMatrix = [propagator]() {
958 return propagator->viewOnPRScalingMatrix();
960 propagatorConnection.getPRScalingCallback = [propagator]() {
961 return propagator->prScalingCallback();
965 // NOLINTNEXTLINE(readability-misleading-indentation)
966 return propagatorConnection;
969 // doxygen is confused by the two definitions
971 template<IntegrationStage integrationStage>
972 ISimulatorElement* Propagator<integrationStage>::getElementPointerImpl(
973 LegacySimulatorData* legacySimulatorData,
974 ModularSimulatorAlgorithmBuilderHelper* builderHelper,
975 StatePropagatorData* statePropagatorData,
976 EnergyData gmx_unused* energyData,
977 FreeEnergyPerturbationData gmx_unused* freeEnergyPerturbationData,
978 GlobalCommunicationHelper gmx_unused* globalCommunicationHelper,
979 const PropagatorTag& propagatorTag,
982 GMX_RELEASE_ASSERT(!(integrationStage == IntegrationStage::ScaleVelocities
983 || integrationStage == IntegrationStage::ScalePositions)
984 || (timestep == 0.0),
985 "Scaling elements don't propagate the system.");
986 auto* element = builderHelper->storeElement(std::make_unique<Propagator<integrationStage>>(
987 timestep, statePropagatorData, legacySimulatorData->mdAtoms, legacySimulatorData->wcycle));
988 auto* propagator = static_cast<Propagator<integrationStage>*>(element);
989 builderHelper->registerPropagator(getConnection<integrationStage>(propagator, propagatorTag));
993 template<IntegrationStage integrationStage>
994 ISimulatorElement* Propagator<integrationStage>::getElementPointerImpl(
995 LegacySimulatorData* legacySimulatorData,
996 ModularSimulatorAlgorithmBuilderHelper* builderHelper,
997 StatePropagatorData* statePropagatorData,
998 EnergyData* energyData,
999 FreeEnergyPerturbationData* freeEnergyPerturbationData,
1000 GlobalCommunicationHelper* globalCommunicationHelper,
1001 const PropagatorTag& propagatorTag)
1004 integrationStage == IntegrationStage::ScaleVelocities
1005 || integrationStage == IntegrationStage::ScalePositions,
1006 "Adding a propagator without time step is only allowed for scaling elements");
1007 return getElementPointerImpl(legacySimulatorData,
1009 statePropagatorData,
1011 freeEnergyPerturbationData,
1012 globalCommunicationHelper,
1018 // Explicit template initializations
1019 template class Propagator<IntegrationStage::PositionsOnly>;
1020 template class Propagator<IntegrationStage::VelocitiesOnly>;
1021 template class Propagator<IntegrationStage::LeapFrog>;
1022 template class Propagator<IntegrationStage::VelocityVerletPositionsAndVelocities>;
1023 template class Propagator<IntegrationStage::ScaleVelocities>;
1024 template class Propagator<IntegrationStage::ScalePositions>;