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