src/gromacs/modularsimulator/propagator.cpp

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