src/gromacs/modularsimulator/propagator.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 propagator element 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 "propagator.h"
  45
  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/mdatom.h"
  53 #include "gromacs/timing/wallcycle.h"
  54 #include "gromacs/utility/fatalerror.h"
  55
  56 #include "statepropagatordata.h"
  57
  58 namespace gmx
  59 {
  60 //! Update velocities
  61 template<NumVelocityScalingValues numVelocityScalingValues, ParrinelloRahmanVelocityScaling parrinelloRahmanVelocityScaling>
  62 static void inline updateVelocities(int         a,
  63                                     real        dt,
  64                                     real        lambda,
  65                                     const rvec* gmx_restrict invMassPerDim,
  66                                     rvec* gmx_restrict v,
  67                                     const rvec* gmx_restrict f,
  68                                     const rvec               diagPR,
  69                                     const matrix             matrixPR)
  70 {
  71     for (int d = 0; d < DIM; d++)
  72     {
  73         // TODO: Extract this into policy classes
  74         if (numVelocityScalingValues != NumVelocityScalingValues::None
  75             && parrinelloRahmanVelocityScaling == ParrinelloRahmanVelocityScaling::No)
  76         {
  77             v[a][d] *= lambda;
  78         }
  79         if (numVelocityScalingValues != NumVelocityScalingValues::None
  80             && parrinelloRahmanVelocityScaling == ParrinelloRahmanVelocityScaling::Diagonal)
  81         {
  82             v[a][d] *= (lambda - diagPR[d]);
  83         }
  84         if (numVelocityScalingValues != NumVelocityScalingValues::None
  85             && parrinelloRahmanVelocityScaling == ParrinelloRahmanVelocityScaling::Full)
  86         {
  87             v[a][d] = lambda * v[a][d] - iprod(matrixPR[d], v[a]);
  88         }
  89         if (numVelocityScalingValues == NumVelocityScalingValues::None
  90             && parrinelloRahmanVelocityScaling == ParrinelloRahmanVelocityScaling::Diagonal)
  91         {
  92             v[a][d] *= (1 - diagPR[d]);
  93         }
  94         if (numVelocityScalingValues == NumVelocityScalingValues::None
  95             && parrinelloRahmanVelocityScaling == ParrinelloRahmanVelocityScaling::Full)
  96         {
  97             v[a][d] -= iprod(matrixPR[d], v[a]);
  98         }
  99         v[a][d] += f[a][d] * invMassPerDim[a][d] * dt;
 100     }
 101 }
 102
 103 //! Update positions
 104 static void inline updatePositions(int         a,
 105                                    real        dt,
 106                                    const rvec* gmx_restrict x,
 107                                    rvec* gmx_restrict xprime,
 108                                    const rvec* gmx_restrict v)
 109 {
 110     for (int d = 0; d < DIM; d++)
 111     {
 112         xprime[a][d] = x[a][d] + v[a][d] * dt;
 113     }
 114 }
 115
 116 //! Propagation (position only)
 117 template<>
 118 template<NumVelocityScalingValues numVelocityScalingValues, ParrinelloRahmanVelocityScaling parrinelloRahmanVelocityScaling>
 119 void Propagator<IntegrationStep::PositionsOnly>::run()
 120 {
 121     wallcycle_start(wcycle_, ewcUPDATE);
 122
 123     auto xp = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
 124     auto x = as_rvec_array(statePropagatorData_->constPreviousPositionsView().paddedArrayRef().data());
 125     auto v = as_rvec_array(statePropagatorData_->constVelocitiesView().paddedArrayRef().data());
 126
 127     int nth    = gmx_omp_nthreads_get(emntUpdate);
 128     int homenr = mdAtoms_->mdatoms()->homenr;
 129
 130 #pragma omp parallel for num_threads(nth) schedule(static) default(none) shared(nth, homenr, x, xp, v)
 131     for (int th = 0; th < nth; th++)
 132     {
 133         try
 134         {
 135             int start_th, end_th;
 136             getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
 137
 138             for (int a = start_th; a < end_th; a++)
 139             {
 140                 updatePositions(a, timestep_, x, xp, v);
 141             }
 142         }
 143         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 144     }
 145     wallcycle_stop(wcycle_, ewcUPDATE);
 146 }
 147
 148 //! Propagation (velocity only)
 149 template<>
 150 template<NumVelocityScalingValues numVelocityScalingValues, ParrinelloRahmanVelocityScaling parrinelloRahmanVelocityScaling>
 151 void Propagator<IntegrationStep::VelocitiesOnly>::run()
 152 {
 153     wallcycle_start(wcycle_, ewcUPDATE);
 154
 155     auto v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
 156     auto f = as_rvec_array(statePropagatorData_->constForcesView().paddedArrayRef().data());
 157     auto invMassPerDim = mdAtoms_->mdatoms()->invMassPerDim;
 158
 159     const real lambda =
 160             (numVelocityScalingValues == NumVelocityScalingValues::Single) ? velocityScaling_[0] : 1.0;
 161
 162     bool doDiagonalScaling = false;
 163     if (parrinelloRahmanVelocityScaling != ParrinelloRahmanVelocityScaling::No)
 164     {
 165         // TODO: Could we know in advance whether the matrix is diagonal?
 166         doDiagonalScaling = (matrixPR[YY][XX] == 0 && matrixPR[ZZ][XX] == 0 && matrixPR[ZZ][YY] == 0);
 167         if (doDiagonalScaling)
 168         {
 169             diagPR[XX] = matrixPR[XX][XX];
 170             diagPR[YY] = matrixPR[YY][YY];
 171             diagPR[ZZ] = matrixPR[ZZ][ZZ];
 172         }
 173     }
 174
 175     int nth    = gmx_omp_nthreads_get(emntUpdate);
 176     int homenr = mdAtoms_->mdatoms()->homenr;
 177
 178 // lambda could be shared, but gcc-8 & gcc-9 don't agree how to write that...
 179 // https://www.gnu.org/software/gcc/gcc-9/porting_to.html -> OpenMP data sharing
 180 #pragma omp parallel for num_threads(nth) schedule(static) default(none) \
 181         shared(nth, homenr, v, f, invMassPerDim, doDiagonalScaling) firstprivate(lambda)
 182     for (int th = 0; th < nth; th++)
 183     {
 184         try
 185         {
 186             int start_th, end_th;
 187             getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
 188
 189             for (int a = start_th; a < end_th; a++)
 190             {
 191                 if (parrinelloRahmanVelocityScaling == ParrinelloRahmanVelocityScaling::No)
 192                 {
 193                     if (numVelocityScalingValues == NumVelocityScalingValues::Multiple)
 194                     {
 195                         updateVelocities<NumVelocityScalingValues::Multiple, ParrinelloRahmanVelocityScaling::No>(
 196                                 a, timestep_, velocityScaling_[mdAtoms_->mdatoms()->cTC[a]],
 197                                 invMassPerDim, v, f, diagPR, matrixPR);
 198                     }
 199                     else
 200                     {
 201                         updateVelocities<numVelocityScalingValues, ParrinelloRahmanVelocityScaling::No>(
 202                                 a, timestep_, lambda, invMassPerDim, v, f, diagPR, matrixPR);
 203                     }
 204                 }
 205                 else
 206                 {
 207                     if (doDiagonalScaling)
 208                     {
 209                         if (numVelocityScalingValues == NumVelocityScalingValues::Multiple)
 210                         {
 211                             updateVelocities<NumVelocityScalingValues::Multiple, ParrinelloRahmanVelocityScaling::Diagonal>(
 212                                     a, timestep_, velocityScaling_[mdAtoms_->mdatoms()->cTC[a]],
 213                                     invMassPerDim, v, f, diagPR, matrixPR);
 214                         }
 215                         else
 216                         {
 217                             updateVelocities<numVelocityScalingValues, ParrinelloRahmanVelocityScaling::Diagonal>(
 218                                     a, timestep_, lambda, invMassPerDim, v, f, diagPR, matrixPR);
 219                         }
 220                     }
 221                     else
 222                     {
 223                         if (numVelocityScalingValues == NumVelocityScalingValues::Multiple)
 224                         {
 225                             updateVelocities<NumVelocityScalingValues::Multiple, ParrinelloRahmanVelocityScaling::Full>(
 226                                     a, timestep_, velocityScaling_[mdAtoms_->mdatoms()->cTC[a]],
 227                                     invMassPerDim, v, f, diagPR, matrixPR);
 228                         }
 229                         else
 230                         {
 231                             updateVelocities<numVelocityScalingValues, ParrinelloRahmanVelocityScaling::Full>(
 232                                     a, timestep_, lambda, invMassPerDim, v, f, diagPR, matrixPR);
 233                         }
 234                     }
 235                 }
 236             }
 237         }
 238         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 239     }
 240     wallcycle_stop(wcycle_, ewcUPDATE);
 241 }
 242
 243 //! Propagation (leapfrog case - position and velocity)
 244 template<>
 245 template<NumVelocityScalingValues numVelocityScalingValues, ParrinelloRahmanVelocityScaling parrinelloRahmanVelocityScaling>
 246 void Propagator<IntegrationStep::LeapFrog>::run()
 247 {
 248     wallcycle_start(wcycle_, ewcUPDATE);
 249
 250     auto xp = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
 251     auto x = as_rvec_array(statePropagatorData_->constPreviousPositionsView().paddedArrayRef().data());
 252     auto v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
 253     auto f = as_rvec_array(statePropagatorData_->constForcesView().paddedArrayRef().data());
 254     auto invMassPerDim = mdAtoms_->mdatoms()->invMassPerDim;
 255
 256     const real lambda =
 257             (numVelocityScalingValues == NumVelocityScalingValues::Single) ? velocityScaling_[0] : 1.0;
 258
 259     bool doDiagonalScaling = false;
 260     if (parrinelloRahmanVelocityScaling != ParrinelloRahmanVelocityScaling::No)
 261     {
 262         // TODO: Could we know in advance whether the matrix is diagonal?
 263         doDiagonalScaling = (matrixPR[YY][XX] == 0 && matrixPR[ZZ][XX] == 0 && matrixPR[ZZ][YY] == 0);
 264         if (doDiagonalScaling)
 265         {
 266             diagPR[XX] = matrixPR[XX][XX];
 267             diagPR[YY] = matrixPR[YY][YY];
 268             diagPR[ZZ] = matrixPR[ZZ][ZZ];
 269         }
 270     }
 271
 272     int nth    = gmx_omp_nthreads_get(emntUpdate);
 273     int homenr = mdAtoms_->mdatoms()->homenr;
 274
 275 // lambda could be shared, but gcc-8 & gcc-9 don't agree how to write that...
 276 // https://www.gnu.org/software/gcc/gcc-9/porting_to.html -> OpenMP data sharing
 277 #pragma omp parallel for num_threads(nth) schedule(static) default(none) \
 278         shared(nth, homenr, x, xp, v, f, invMassPerDim, doDiagonalScaling) firstprivate(lambda)
 279     for (int th = 0; th < nth; th++)
 280     {
 281         try
 282         {
 283             int start_th, end_th;
 284             getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
 285
 286             for (int a = start_th; a < end_th; a++)
 287             {
 288                 if (parrinelloRahmanVelocityScaling == ParrinelloRahmanVelocityScaling::No)
 289                 {
 290                     if (numVelocityScalingValues == NumVelocityScalingValues::Multiple)
 291                     {
 292                         updateVelocities<NumVelocityScalingValues::Multiple, ParrinelloRahmanVelocityScaling::No>(
 293                                 a, timestep_, velocityScaling_[mdAtoms_->mdatoms()->cTC[a]],
 294                                 invMassPerDim, v, f, diagPR, matrixPR);
 295                     }
 296                     else
 297                     {
 298                         updateVelocities<numVelocityScalingValues, ParrinelloRahmanVelocityScaling::No>(
 299                                 a, timestep_, lambda, invMassPerDim, v, f, diagPR, matrixPR);
 300                     }
 301                 }
 302                 else
 303                 {
 304                     if (doDiagonalScaling)
 305                     {
 306                         if (numVelocityScalingValues == NumVelocityScalingValues::Multiple)
 307                         {
 308                             updateVelocities<NumVelocityScalingValues::Multiple, ParrinelloRahmanVelocityScaling::Diagonal>(
 309                                     a, timestep_, velocityScaling_[mdAtoms_->mdatoms()->cTC[a]],
 310                                     invMassPerDim, v, f, diagPR, matrixPR);
 311                         }
 312                         else
 313                         {
 314                             updateVelocities<numVelocityScalingValues, ParrinelloRahmanVelocityScaling::Diagonal>(
 315                                     a, timestep_, lambda, invMassPerDim, v, f, diagPR, matrixPR);
 316                         }
 317                     }
 318                     else
 319                     {
 320                         if (numVelocityScalingValues == NumVelocityScalingValues::Multiple)
 321                         {
 322                             updateVelocities<NumVelocityScalingValues::Multiple, ParrinelloRahmanVelocityScaling::Full>(
 323                                     a, timestep_, velocityScaling_[mdAtoms_->mdatoms()->cTC[a]],
 324                                     invMassPerDim, v, f, diagPR, matrixPR);
 325                         }
 326                         else
 327                         {
 328                             updateVelocities<numVelocityScalingValues, ParrinelloRahmanVelocityScaling::Full>(
 329                                     a, timestep_, lambda, invMassPerDim, v, f, diagPR, matrixPR);
 330                         }
 331                     }
 332                 }
 333                 updatePositions(a, timestep_, x, xp, v);
 334             }
 335         }
 336         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 337     }
 338     wallcycle_stop(wcycle_, ewcUPDATE);
 339 }
 340
 341 //! Propagation (velocity verlet stage 2 - velocity and position)
 342 template<>
 343 template<NumVelocityScalingValues numVelocityScalingValues, ParrinelloRahmanVelocityScaling parrinelloRahmanVelocityScaling>
 344 void Propagator<IntegrationStep::VelocityVerletPositionsAndVelocities>::run()
 345 {
 346     wallcycle_start(wcycle_, ewcUPDATE);
 347
 348     auto xp = as_rvec_array(statePropagatorData_->positionsView().paddedArrayRef().data());
 349     auto x = as_rvec_array(statePropagatorData_->constPreviousPositionsView().paddedArrayRef().data());
 350     auto v = as_rvec_array(statePropagatorData_->velocitiesView().paddedArrayRef().data());
 351     auto f = as_rvec_array(statePropagatorData_->constForcesView().paddedArrayRef().data());
 352     auto invMassPerDim = mdAtoms_->mdatoms()->invMassPerDim;
 353
 354     int nth    = gmx_omp_nthreads_get(emntUpdate);
 355     int homenr = mdAtoms_->mdatoms()->homenr;
 356
 357     const real lambda =
 358             (numVelocityScalingValues == NumVelocityScalingValues::Single) ? velocityScaling_[0] : 1.0;
 359
 360     bool doDiagonalScaling = false;
 361     if (parrinelloRahmanVelocityScaling != ParrinelloRahmanVelocityScaling::No)
 362     {
 363         // TODO: Could we know in advance whether the matrix is diagonal?
 364         doDiagonalScaling = (matrixPR[YY][XX] == 0 && matrixPR[ZZ][XX] == 0 && matrixPR[ZZ][YY] == 0);
 365         if (doDiagonalScaling)
 366         {
 367             diagPR[XX] = matrixPR[XX][XX];
 368             diagPR[YY] = matrixPR[YY][YY];
 369             diagPR[ZZ] = matrixPR[ZZ][ZZ];
 370         }
 371     }
 372
 373 // lambda could be shared, but gcc-8 & gcc-9 don't agree how to write that...
 374 // https://www.gnu.org/software/gcc/gcc-9/porting_to.html -> OpenMP data sharing
 375 #pragma omp parallel for num_threads(nth) schedule(static) default(none) \
 376         shared(nth, homenr, x, xp, v, f, invMassPerDim, doDiagonalScaling) firstprivate(lambda)
 377     for (int th = 0; th < nth; th++)
 378     {
 379         try
 380         {
 381             int start_th, end_th;
 382             getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
 383
 384             for (int a = start_th; a < end_th; a++)
 385             {
 386                 if (parrinelloRahmanVelocityScaling == ParrinelloRahmanVelocityScaling::No)
 387                 {
 388                     if (numVelocityScalingValues == NumVelocityScalingValues::Multiple)
 389                     {
 390                         updateVelocities<NumVelocityScalingValues::Multiple, ParrinelloRahmanVelocityScaling::No>(
 391                                 a, 0.5 * timestep_, velocityScaling_[mdAtoms_->mdatoms()->cTC[a]],
 392                                 invMassPerDim, v, f, diagPR, matrixPR);
 393                     }
 394                     else
 395                     {
 396                         updateVelocities<numVelocityScalingValues, ParrinelloRahmanVelocityScaling::No>(
 397                                 a, 0.5 * timestep_, lambda, invMassPerDim, v, f, diagPR, matrixPR);
 398                     }
 399                 }
 400                 else
 401                 {
 402                     if (doDiagonalScaling)
 403                     {
 404                         if (numVelocityScalingValues == NumVelocityScalingValues::Multiple)
 405                         {
 406                             updateVelocities<NumVelocityScalingValues::Multiple, ParrinelloRahmanVelocityScaling::Diagonal>(
 407                                     a, 0.5 * timestep_, velocityScaling_[mdAtoms_->mdatoms()->cTC[a]],
 408                                     invMassPerDim, v, f, diagPR, matrixPR);
 409                         }
 410                         else
 411                         {
 412                             updateVelocities<numVelocityScalingValues, ParrinelloRahmanVelocityScaling::Diagonal>(
 413                                     a, 0.5 * timestep_, lambda, invMassPerDim, v, f, diagPR, matrixPR);
 414                         }
 415                     }
 416                     else
 417                     {
 418                         if (numVelocityScalingValues == NumVelocityScalingValues::Multiple)
 419                         {
 420                             updateVelocities<NumVelocityScalingValues::Multiple, ParrinelloRahmanVelocityScaling::Full>(
 421                                     a, 0.5 * timestep_, velocityScaling_[mdAtoms_->mdatoms()->cTC[a]],
 422                                     invMassPerDim, v, f, diagPR, matrixPR);
 423                         }
 424                         else
 425                         {
 426                             updateVelocities<numVelocityScalingValues, ParrinelloRahmanVelocityScaling::Full>(
 427                                     a, 0.5 * timestep_, lambda, invMassPerDim, v, f, diagPR, matrixPR);
 428                         }
 429                     }
 430                 }
 431                 updatePositions(a, timestep_, x, xp, v);
 432             }
 433         }
 434         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
 435     }
 436     wallcycle_stop(wcycle_, ewcUPDATE);
 437 }
 438
 439 template<IntegrationStep algorithm>
 440 Propagator<algorithm>::Propagator(double               timestep,
 441                                   StatePropagatorData* statePropagatorData,
 442                                   const MDAtoms*       mdAtoms,
 443                                   gmx_wallcycle*       wcycle) :
 444     timestep_(timestep),
 445     statePropagatorData_(statePropagatorData),
 446     doSingleVelocityScaling(false),
 447     doGroupVelocityScaling(false),
 448     scalingStepVelocity_(-1),
 449     scalingStepPR_(-1),
 450     mdAtoms_(mdAtoms),
 451     wcycle_(wcycle)
 452 {
 453     clear_rvec(diagPR);
 454     clear_mat(matrixPR);
 455 }
 456
 457 template<IntegrationStep algorithm>
 458 void Propagator<algorithm>::scheduleTask(Step gmx_unused step,
 459                                          Time gmx_unused               time,
 460                                          const RegisterRunFunctionPtr& registerRunFunction)
 461 {
 462     const bool doSingleVScalingThisStep = (doSingleVelocityScaling && (step == scalingStepVelocity_));
 463     const bool doGroupVScalingThisStep = (doGroupVelocityScaling && (step == scalingStepVelocity_));
 464
 465     const bool doParrinelloRahmanThisStep = (step == scalingStepPR_);
 466
 467     if (doSingleVScalingThisStep)
 468     {
 469         if (doParrinelloRahmanThisStep)
 470         {
 471             (*registerRunFunction)(std::make_unique<SimulatorRunFunction>([this]() {
 472                 run<NumVelocityScalingValues::Single, ParrinelloRahmanVelocityScaling::Full>();
 473             }));
 474         }
 475         else
 476         {
 477             (*registerRunFunction)(std::make_unique<SimulatorRunFunction>([this]() {
 478                 run<NumVelocityScalingValues::Single, ParrinelloRahmanVelocityScaling::No>();
 479             }));
 480         }
 481     }
 482     else if (doGroupVScalingThisStep)
 483     {
 484         if (doParrinelloRahmanThisStep)
 485         {
 486             (*registerRunFunction)(std::make_unique<SimulatorRunFunction>([this]() {
 487                 run<NumVelocityScalingValues::Multiple, ParrinelloRahmanVelocityScaling::Full>();
 488             }));
 489         }
 490         else
 491         {
 492             (*registerRunFunction)(std::make_unique<SimulatorRunFunction>([this]() {
 493                 run<NumVelocityScalingValues::Multiple, ParrinelloRahmanVelocityScaling::No>();
 494             }));
 495         }
 496     }
 497     else
 498     {
 499         if (doParrinelloRahmanThisStep)
 500         {
 501             (*registerRunFunction)(std::make_unique<SimulatorRunFunction>([this]() {
 502                 run<NumVelocityScalingValues::None, ParrinelloRahmanVelocityScaling::Full>();
 503             }));
 504         }
 505         else
 506         {
 507             (*registerRunFunction)(std::make_unique<SimulatorRunFunction>([this]() {
 508                 run<NumVelocityScalingValues::None, ParrinelloRahmanVelocityScaling::No>();
 509             }));
 510         }
 511     }
 512 }
 513
 514 template<IntegrationStep algorithm>
 515 void Propagator<algorithm>::setNumVelocityScalingVariables(int numVelocityScalingVariables)
 516 {
 517     if (algorithm == IntegrationStep::PositionsOnly)
 518     {
 519         gmx_fatal(FARGS, "Velocity scaling not implemented for IntegrationStep::PositionsOnly.");
 520     }
 521     GMX_ASSERT(velocityScaling_.empty(),
 522                "Number of velocity scaling variables cannot be changed once set.");
 523
 524     velocityScaling_.resize(numVelocityScalingVariables, 1.);
 525     doSingleVelocityScaling = numVelocityScalingVariables == 1;
 526     doGroupVelocityScaling  = numVelocityScalingVariables > 1;
 527 }
 528
 529 template<IntegrationStep algorithm>
 530 ArrayRef<real> Propagator<algorithm>::viewOnVelocityScaling()
 531 {
 532     if (algorithm == IntegrationStep::PositionsOnly)
 533     {
 534         gmx_fatal(FARGS, "Velocity scaling not implemented for IntegrationStep::PositionsOnly.");
 535     }
 536     GMX_ASSERT(!velocityScaling_.empty(), "Number of velocity scaling variables not set.");
 537
 538     return velocityScaling_;
 539 }
 540
 541 template<IntegrationStep algorithm>
 542 std::unique_ptr<std::function<void(Step)>> Propagator<algorithm>::velocityScalingCallback()
 543 {
 544     if (algorithm == IntegrationStep::PositionsOnly)
 545     {
 546         gmx_fatal(FARGS, "Velocity scaling not implemented for IntegrationStep::PositionsOnly.");
 547     }
 548
 549     return std::make_unique<PropagatorCallback>([this](Step step) { scalingStepVelocity_ = step; });
 550 }
 551
 552 template<IntegrationStep algorithm>
 553 ArrayRef<rvec> Propagator<algorithm>::viewOnPRScalingMatrix()
 554 {
 555     GMX_RELEASE_ASSERT(
 556             algorithm != IntegrationStep::PositionsOnly,
 557             "Parrinello-Rahman scaling not implemented for IntegrationStep::PositionsOnly.");
 558
 559     clear_mat(matrixPR);
 560     // gcc-5 needs this to be explicit (all other tested compilers would be ok
 561     // with simply returning matrixPR)
 562     return ArrayRef<rvec>(matrixPR);
 563 }
 564
 565 template<IntegrationStep algorithm>
 566 PropagatorCallbackPtr Propagator<algorithm>::prScalingCallback()
 567 {
 568     GMX_RELEASE_ASSERT(
 569             algorithm != IntegrationStep::PositionsOnly,
 570             "Parrinello-Rahman scaling not implemented for IntegrationStep::PositionsOnly.");
 571
 572     return std::make_unique<PropagatorCallback>([this](Step step) { scalingStepPR_ = step; });
 573 }
 574
 575 //! Explicit template initialization
 576 //! @{
 577 template class Propagator<IntegrationStep::PositionsOnly>;
 578 template class Propagator<IntegrationStep::VelocitiesOnly>;
 579 template class Propagator<IntegrationStep::LeapFrog>;
 580 template class Propagator<IntegrationStep::VelocityVerletPositionsAndVelocities>;
 581 //! @}
 582
 583 } // namespace gmx