src/gromacs/modularsimulator/propagator.cpp

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