src/gromacs/mdlib/update.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
   5  * Copyright (c) 2001-2004, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 #include "gmxpre.h"
  38
  39 #include "update.h"
  40
  41 #include <cmath>
  42 #include <cstdio>
  43
  44 #include <algorithm>
  45 #include <memory>
  46
  47 #include "gromacs/domdec/domdec_struct.h"
  48 #include "gromacs/fileio/confio.h"
  49 #include "gromacs/gmxlib/network.h"
  50 #include "gromacs/gmxlib/nrnb.h"
  51 #include "gromacs/listed_forces/disre.h"
  52 #include "gromacs/listed_forces/orires.h"
  53 #include "gromacs/math/functions.h"
  54 #include "gromacs/math/invertmatrix.h"
  55 #include "gromacs/math/paddedvector.h"
  56 #include "gromacs/math/units.h"
  57 #include "gromacs/math/vec.h"
  58 #include "gromacs/math/vecdump.h"
  59 #include "gromacs/mdlib/boxdeformation.h"
  60 #include "gromacs/mdlib/constr.h"
  61 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  62 #include "gromacs/mdlib/mdatoms.h"
  63 #include "gromacs/mdlib/stat.h"
  64 #include "gromacs/mdlib/tgroup.h"
  65 #include "gromacs/mdtypes/commrec.h"
  66 #include "gromacs/mdtypes/group.h"
  67 #include "gromacs/mdtypes/inputrec.h"
  68 #include "gromacs/mdtypes/md_enums.h"
  69 #include "gromacs/mdtypes/state.h"
  70 #include "gromacs/pbcutil/boxutilities.h"
  71 #include "gromacs/pbcutil/mshift.h"
  72 #include "gromacs/pbcutil/pbc.h"
  73 #include "gromacs/pulling/pull.h"
  74 #include "gromacs/random/tabulatednormaldistribution.h"
  75 #include "gromacs/random/threefry.h"
  76 #include "gromacs/simd/simd.h"
  77 #include "gromacs/timing/wallcycle.h"
  78 #include "gromacs/topology/atoms.h"
  79 #include "gromacs/utility/exceptions.h"
  80 #include "gromacs/utility/fatalerror.h"
  81 #include "gromacs/utility/futil.h"
  82 #include "gromacs/utility/gmxassert.h"
  83 #include "gromacs/utility/gmxomp.h"
  84 #include "gromacs/utility/smalloc.h"
  85
  86 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
  87
  88 struct gmx_sd_const_t
  89 {
  90     double em = 0;
  91 };
  92
  93 struct gmx_sd_sigma_t
  94 {
  95     real V = 0;
  96 };
  97
  98 struct gmx_stochd_t
  99 {
 100     /* BD stuff */
 101     std::vector<real>           bd_rf;
 102     /* SD stuff */
 103     std::vector<gmx_sd_const_t> sdc;
 104     std::vector<gmx_sd_sigma_t> sdsig;
 105     /* andersen temperature control stuff */
 106     std::vector<bool>           randomize_group;
 107     std::vector<real>           boltzfac;
 108
 109     explicit gmx_stochd_t(const t_inputrec *ir);
 110 };
 111
 112 //! pImpled implementation for Update
 113 class Update::Impl
 114 {
 115     public:
 116         //! Constructor
 117         Impl(const t_inputrec    *ir, BoxDeformation *boxDeformation);
 118         //! Destructor
 119         ~Impl() = default;
 120         //! stochastic dynamics struct
 121         std::unique_ptr<gmx_stochd_t> sd;
 122         //! xprime for constraint algorithms
 123         PaddedVector<RVec>            xp;
 124         //! Box deformation handler (or nullptr if inactive).
 125         BoxDeformation               *deform = nullptr;
 126 };
 127
 128 Update::Update(const t_inputrec    *ir, BoxDeformation *boxDeformation)
 129     : impl_(new Impl(ir, boxDeformation))
 130 {};
 131
 132 Update::~Update()
 133 {};
 134
 135 gmx_stochd_t* Update::sd() const
 136 {
 137     return impl_->sd.get();
 138 }
 139
 140 PaddedVector<RVec> * Update::xp()
 141 {
 142     return &impl_->xp;
 143 }
 144
 145 BoxDeformation * Update::deform() const
 146 {
 147     return impl_->deform;
 148 }
 149
 150 /*! \brief Sets the velocities of virtual sites to zero */
 151 static void clearVsiteVelocities(int                   start,
 152                                  int                   nrend,
 153                                  const unsigned short *particleType,
 154                                  rvec * gmx_restrict   v)
 155 {
 156     for (int a = start; a < nrend; a++)
 157     {
 158         if (particleType[a] == eptVSite)
 159         {
 160             clear_rvec(v[a]);
 161         }
 162     }
 163 }
 164
 165 /*! \brief Sets the number of different temperature coupling values */
 166 enum class NumTempScaleValues
 167 {
 168     single,   //!< Single T-scaling value (either one group or all values =1)
 169     multiple  //!< Multiple T-scaling values, need to use T-group indices
 170 };
 171
 172 /*! \brief Sets if to apply no or diagonal Parrinello-Rahman pressure scaling
 173  *
 174  * Note that this enum is only used in updateMDLeapfrogSimple(), which does
 175  * not handle fully anistropic Parrinello-Rahman scaling, so we only have
 176  * options \p no and \p diagonal here and no anistropic option.
 177  */
 178 enum class ApplyParrinelloRahmanVScaling
 179 {
 180     no,       //!< Do not apply velocity scaling (not a PR-coupling run or step)
 181     diagonal  //!< Apply velocity scaling using a diagonal matrix
 182 };
 183
 184 /*! \brief Integrate using leap-frog with T-scaling and optionally diagonal Parrinello-Rahman p-coupling
 185  *
 186  * \tparam       numTempScaleValues     The number of different T-couple values
 187  * \tparam       applyPRVScaling        Apply Parrinello-Rahman velocity scaling
 188  * \param[in]    start                  Index of first atom to update
 189  * \param[in]    nrend                  Last atom to update: \p nrend - 1
 190  * \param[in]    dt                     The time step
 191  * \param[in]    dtPressureCouple       Time step for pressure coupling
 192  * \param[in]    invMassPerDim          1/mass per atom and dimension
 193  * \param[in]    tcstat                 Temperature coupling information
 194  * \param[in]    cTC                    T-coupling group index per atom
 195  * \param[in]    pRVScaleMatrixDiagonal Parrinello-Rahman v-scale matrix diagonal
 196  * \param[in]    x                      Input coordinates
 197  * \param[out]   xprime                 Updated coordinates
 198  * \param[inout] v                      Velocities
 199  * \param[in]    f                      Forces
 200  *
 201  * We expect this template to get good SIMD acceleration by most compilers,
 202  * unlike the more complex general template.
 203  * Note that we might get even better SIMD acceleration when we introduce
 204  * aligned (and padded) memory, possibly with some hints for the compilers.
 205  */
 206 template<NumTempScaleValues            numTempScaleValues,
 207          ApplyParrinelloRahmanVScaling applyPRVScaling>
 208 static void
 209 updateMDLeapfrogSimple(int                               start,
 210                        int                               nrend,
 211                        real                              dt,
 212                        real                              dtPressureCouple,
 213                        const rvec * gmx_restrict         invMassPerDim,
 214                        gmx::ArrayRef<const t_grp_tcstat> tcstat,
 215                        const unsigned short            * cTC,
 216                        const rvec                        pRVScaleMatrixDiagonal,
 217                        const rvec * gmx_restrict         x,
 218                        rvec       * gmx_restrict         xprime,
 219                        rvec       * gmx_restrict         v,
 220                        const rvec * gmx_restrict         f)
 221 {
 222     real lambdaGroup;
 223
 224     if (numTempScaleValues == NumTempScaleValues::single)
 225     {
 226         lambdaGroup = tcstat[0].lambda;
 227     }
 228
 229     for (int a = start; a < nrend; a++)
 230     {
 231         if (numTempScaleValues == NumTempScaleValues::multiple)
 232         {
 233             lambdaGroup = tcstat[cTC[a]].lambda;
 234         }
 235
 236         for (int d = 0; d < DIM; d++)
 237         {
 238             /* Note that using rvec invMassPerDim results in more efficient
 239              * SIMD code, but this increases the cache pressure.
 240              * For large systems with PME on the CPU this slows down the
 241              * (then already slow) update by 20%. If all data remains in cache,
 242              * using rvec is much faster.
 243              */
 244             real vNew = lambdaGroup*v[a][d] + f[a][d]*invMassPerDim[a][d]*dt;
 245
 246             if (applyPRVScaling == ApplyParrinelloRahmanVScaling::diagonal)
 247             {
 248                 vNew -= dtPressureCouple*pRVScaleMatrixDiagonal[d]*v[a][d];
 249             }
 250             v[a][d]      = vNew;
 251             xprime[a][d] = x[a][d] + vNew*dt;
 252         }
 253     }
 254 }
 255
 256 #if GMX_SIMD && GMX_SIMD_HAVE_REAL
 257 #define GMX_HAVE_SIMD_UPDATE 1
 258 #else
 259 #define GMX_HAVE_SIMD_UPDATE 0
 260 #endif
 261
 262 #if GMX_HAVE_SIMD_UPDATE
 263
 264 /*! \brief Load (aligned) the contents of GMX_SIMD_REAL_WIDTH rvec elements sequentially into 3 SIMD registers
 265  *
 266  * The loaded output is:
 267  * \p r0: { r[index][XX], r[index][YY], ... }
 268  * \p r1: { ... }
 269  * \p r2: { ..., r[index+GMX_SIMD_REAL_WIDTH-1][YY], r[index+GMX_SIMD_REAL_WIDTH-1][ZZ] }
 270  *
 271  * \param[in]  r      Real to an rvec array, has to be aligned to SIMD register width
 272  * \param[in]  index  Index of the first rvec triplet of reals to load
 273  * \param[out] r0     Pointer to first SIMD register
 274  * \param[out] r1     Pointer to second SIMD register
 275  * \param[out] r2     Pointer to third SIMD register
 276  */
 277 static inline void simdLoadRvecs(const rvec *r,
 278                                  int         index,
 279                                  SimdReal   *r0,
 280                                  SimdReal   *r1,
 281                                  SimdReal   *r2)
 282 {
 283     const real *realPtr = r[index];
 284
 285     GMX_ASSERT(isSimdAligned(realPtr), "Pointer should be SIMD aligned");
 286
 287     *r0 = simdLoad(realPtr + 0*GMX_SIMD_REAL_WIDTH);
 288     *r1 = simdLoad(realPtr + 1*GMX_SIMD_REAL_WIDTH);
 289     *r2 = simdLoad(realPtr + 2*GMX_SIMD_REAL_WIDTH);
 290 }
 291
 292 /*! \brief Store (aligned) 3 SIMD registers sequentially to GMX_SIMD_REAL_WIDTH rvec elements
 293  *
 294  * The stored output is:
 295  * \p r[index] = { { r0[0], r0[1], ... }
 296  * ...
 297  * \p r[index+GMX_SIMD_REAL_WIDTH-1] =  { ... , r2[GMX_SIMD_REAL_WIDTH-2], r2[GMX_SIMD_REAL_WIDTH-1] }
 298  *
 299  * \param[out] r      Pointer to an rvec array, has to be aligned to SIMD register width
 300  * \param[in]  index  Index of the first rvec triplet of reals to store to
 301  * \param[in]  r0     First SIMD register
 302  * \param[in]  r1     Second SIMD register
 303  * \param[in]  r2     Third SIMD register
 304  */
 305 static inline void simdStoreRvecs(rvec     *r,
 306                                   int       index,
 307                                   SimdReal  r0,
 308                                   SimdReal  r1,
 309                                   SimdReal  r2)
 310 {
 311     real *realPtr = r[index];
 312
 313     GMX_ASSERT(isSimdAligned(realPtr), "Pointer should be SIMD aligned");
 314
 315     store(realPtr + 0*GMX_SIMD_REAL_WIDTH, r0);
 316     store(realPtr + 1*GMX_SIMD_REAL_WIDTH, r1);
 317     store(realPtr + 2*GMX_SIMD_REAL_WIDTH, r2);
 318 }
 319
 320 /*! \brief Integrate using leap-frog with single group T-scaling and SIMD
 321  *
 322  * \param[in]    start                  Index of first atom to update
 323  * \param[in]    nrend                  Last atom to update: \p nrend - 1
 324  * \param[in]    dt                     The time step
 325  * \param[in]    invMass                1/mass per atom
 326  * \param[in]    tcstat                 Temperature coupling information
 327  * \param[in]    x                      Input coordinates
 328  * \param[out]   xprime                 Updated coordinates
 329  * \param[inout] v                      Velocities
 330  * \param[in]    f                      Forces
 331  */
 332 static void
 333 updateMDLeapfrogSimpleSimd(int                               start,
 334                            int                               nrend,
 335                            real                              dt,
 336                            const real * gmx_restrict         invMass,
 337                            gmx::ArrayRef<const t_grp_tcstat> tcstat,
 338                            const rvec * gmx_restrict         x,
 339                            rvec       * gmx_restrict         xprime,
 340                            rvec       * gmx_restrict         v,
 341                            const rvec * gmx_restrict         f)
 342 {
 343     SimdReal                  timestep(dt);
 344     SimdReal                  lambdaSystem(tcstat[0].lambda);
 345
 346     /* We declare variables here, since code is often slower when declaring them inside the loop */
 347
 348     /* Note: We should implement a proper PaddedVector, so we don't need this check */
 349     GMX_ASSERT(isSimdAligned(invMass), "invMass should be aligned");
 350
 351     for (int a = start; a < nrend; a += GMX_SIMD_REAL_WIDTH)
 352     {
 353         SimdReal invMass0, invMass1, invMass2;
 354         expandScalarsToTriplets(simdLoad(invMass + a),
 355                                 &invMass0, &invMass1, &invMass2);
 356
 357         SimdReal v0, v1, v2;
 358         SimdReal f0, f1, f2;
 359         simdLoadRvecs(v, a, &v0, &v1, &v2);
 360         simdLoadRvecs(f, a, &f0, &f1, &f2);
 361
 362         v0 = fma(f0*invMass0, timestep, lambdaSystem*v0);
 363         v1 = fma(f1*invMass1, timestep, lambdaSystem*v1);
 364         v2 = fma(f2*invMass2, timestep, lambdaSystem*v2);
 365
 366         simdStoreRvecs(v, a, v0, v1, v2);
 367
 368         SimdReal x0, x1, x2;
 369         simdLoadRvecs(x, a, &x0, &x1, &x2);
 370
 371         SimdReal xprime0 = fma(v0, timestep, x0);
 372         SimdReal xprime1 = fma(v1, timestep, x1);
 373         SimdReal xprime2 = fma(v2, timestep, x2);
 374
 375         simdStoreRvecs(xprime, a, xprime0, xprime1, xprime2);
 376     }
 377 }
 378
 379 #endif // GMX_HAVE_SIMD_UPDATE
 380
 381 /*! \brief Sets the NEMD acceleration type */
 382 enum class AccelerationType
 383 {
 384     none, group, cosine
 385 };
 386
 387 /*! \brief Integrate using leap-frog with support for everything
 388  *
 389  * \tparam       accelerationType       Type of NEMD acceleration
 390  * \param[in]    start                  Index of first atom to update
 391  * \param[in]    nrend                  Last atom to update: \p nrend - 1
 392  * \param[in]    doNoseHoover           If to apply Nose-Hoover T-coupling
 393  * \param[in]    dt                     The time step
 394  * \param[in]    dtPressureCouple       Time step for pressure coupling, is 0 when no pressure coupling should be applied at this step
 395  * \param[in]    ir                     The input parameter record
 396  * \param[in]    md                     Atom properties
 397  * \param[in]    ekind                  Kinetic energy data
 398  * \param[in]    box                    The box dimensions
 399  * \param[in]    x                      Input coordinates
 400  * \param[out]   xprime                 Updated coordinates
 401  * \param[inout] v                      Velocities
 402  * \param[in]    f                      Forces
 403  * \param[in]    nh_vxi                 Nose-Hoover velocity scaling factors
 404  * \param[in]    M                      Parrinello-Rahman scaling matrix
 405  */
 406 template<AccelerationType accelerationType>
 407 static void
 408 updateMDLeapfrogGeneral(int                         start,
 409                         int                         nrend,
 410                         bool                        doNoseHoover,
 411                         real                        dt,
 412                         real                        dtPressureCouple,
 413                         const t_inputrec          * ir,
 414                         const t_mdatoms           * md,
 415                         const gmx_ekindata_t      * ekind,
 416                         const matrix                box,
 417                         const rvec   * gmx_restrict x,
 418                         rvec         * gmx_restrict xprime,
 419                         rvec         * gmx_restrict v,
 420                         const rvec   * gmx_restrict f,
 421                         const double * gmx_restrict nh_vxi,
 422                         const matrix                M)
 423 {
 424     /* This is a version of the leap-frog integrator that supports
 425      * all combinations of T-coupling, P-coupling and NEMD.
 426      * Nose-Hoover uses the reversible leap-frog integrator from
 427      * Holian et al. Phys Rev E 52(3) : 2338, 1995
 428      */
 429
 430     gmx::ArrayRef<const t_grp_tcstat> tcstat        = ekind->tcstat;
 431     gmx::ArrayRef<const t_grp_acc>    grpstat       = ekind->grpstat;
 432     const unsigned short            * cTC           = md->cTC;
 433     const unsigned short            * cACC          = md->cACC;
 434     const rvec                      * accel         = ir->opts.acc;
 435
 436     const rvec * gmx_restrict         invMassPerDim = md->invMassPerDim;
 437
 438     /* Initialize group values, changed later when multiple groups are used */
 439     int  ga       = 0;
 440     int  gt       = 0;
 441     real factorNH = 0;
 442
 443     real omega_Z  = 2*static_cast<real>(M_PI)/box[ZZ][ZZ];
 444
 445     for (int n = start; n < nrend; n++)
 446     {
 447         if (cTC)
 448         {
 449             gt  = cTC[n];
 450         }
 451         real lg = tcstat[gt].lambda;
 452
 453         rvec vRel;
 454         real cosineZ, vCosine;
 455 #ifdef __INTEL_COMPILER
 456 #pragma warning( disable : 280 )
 457 #endif
 458         switch (accelerationType)
 459         {
 460             case AccelerationType::none:
 461                 copy_rvec(v[n], vRel);
 462                 break;
 463             case AccelerationType::group:
 464                 if (cACC)
 465                 {
 466                     ga = cACC[n];
 467                 }
 468                 /* Avoid scaling the group velocity */
 469                 rvec_sub(v[n], grpstat[ga].u, vRel);
 470                 break;
 471             case AccelerationType::cosine:
 472                 cosineZ = std::cos(x[n][ZZ]*omega_Z);
 473                 vCosine = cosineZ*ekind->cosacc.vcos;
 474                 /* Avoid scaling the cosine profile velocity */
 475                 copy_rvec(v[n], vRel);
 476                 vRel[XX] -= vCosine;
 477                 break;
 478         }
 479
 480         if (doNoseHoover)
 481         {
 482             /* Here we account for multiple time stepping, by increasing
 483              * the Nose-Hoover correction by nsttcouple
 484              */
 485             factorNH = 0.5*ir->nsttcouple*dt*nh_vxi[gt];
 486         }
 487
 488         for (int d = 0; d < DIM; d++)
 489         {
 490             real vNew =
 491                 (lg*vRel[d] + (f[n][d]*invMassPerDim[n][d]*dt - factorNH*vRel[d]
 492                                - dtPressureCouple*iprod(M[d], vRel)))/(1 + factorNH);
 493             switch (accelerationType)
 494             {
 495                 case AccelerationType::none:
 496                     break;
 497                 case AccelerationType::group:
 498                     /* Add back the mean velocity and apply acceleration */
 499                     vNew += grpstat[ga].u[d] + accel[ga][d]*dt;
 500                     break;
 501                 case AccelerationType::cosine:
 502                     if (d == XX)
 503                     {
 504                         /* Add back the mean velocity and apply acceleration */
 505                         vNew += vCosine + cosineZ*ekind->cosacc.cos_accel*dt;
 506                     }
 507                     break;
 508             }
 509             v[n][d]       = vNew;
 510             xprime[n][d]  = x[n][d] + vNew*dt;
 511         }
 512     }
 513 }
 514
 515 /*! \brief Handles the Leap-frog MD x and v integration */
 516 static void do_update_md(int                         start,
 517                          int                         nrend,
 518                          int64_t                     step,
 519                          real                        dt,
 520                          const t_inputrec          * ir,
 521                          const t_mdatoms           * md,
 522                          const gmx_ekindata_t      * ekind,
 523                          const matrix                box,
 524                          const rvec   * gmx_restrict x,
 525                          rvec         * gmx_restrict xprime,
 526                          rvec         * gmx_restrict v,
 527                          const rvec   * gmx_restrict f,
 528                          const double * gmx_restrict nh_vxi,
 529                          const matrix                M)
 530 {
 531     GMX_ASSERT(nrend == start || xprime != x, "For SIMD optimization certain compilers need to have xprime != x");
 532     GMX_ASSERT(ir->eI == eiMD, "Leap-frog integrator was called while another integrator was requested");
 533
 534     /* Note: Berendsen pressure scaling is handled after do_update_md() */
 535     bool doTempCouple       = (ir->etc != etcNO && do_per_step(step + ir->nsttcouple - 1, ir->nsttcouple));
 536     bool doNoseHoover       = (ir->etc == etcNOSEHOOVER && doTempCouple);
 537     bool doParrinelloRahman = (ir->epc == epcPARRINELLORAHMAN && do_per_step(step + ir->nstpcouple - 1, ir->nstpcouple));
 538     bool doPROffDiagonal    = (doParrinelloRahman && (M[YY][XX] != 0 || M[ZZ][XX] != 0 || M[ZZ][YY] != 0));
 539
 540     real dtPressureCouple   = (doParrinelloRahman ? ir->nstpcouple*dt : 0);
 541
 542     /* NEMD (also cosine) acceleration is applied in updateMDLeapFrogGeneral */
 543     bool doAcceleration     = (ekind->bNEMD || ekind->cosacc.cos_accel != 0);
 544
 545     if (doNoseHoover || doPROffDiagonal || doAcceleration)
 546     {
 547         matrix stepM;
 548         if (!doParrinelloRahman)
 549         {
 550             /* We should not apply PR scaling at this step */
 551             clear_mat(stepM);
 552         }
 553         else
 554         {
 555             copy_mat(M, stepM);
 556         }
 557
 558         if (!doAcceleration)
 559         {
 560             updateMDLeapfrogGeneral<AccelerationType::none>
 561                 (start, nrend, doNoseHoover, dt, dtPressureCouple,
 562                 ir, md, ekind, box, x, xprime, v, f, nh_vxi, stepM);
 563         }
 564         else if (ekind->bNEMD)
 565         {
 566             updateMDLeapfrogGeneral<AccelerationType::group>
 567                 (start, nrend, doNoseHoover, dt, dtPressureCouple,
 568                 ir, md, ekind, box, x, xprime, v, f, nh_vxi, stepM);
 569         }
 570         else
 571         {
 572             updateMDLeapfrogGeneral<AccelerationType::cosine>
 573                 (start, nrend, doNoseHoover, dt, dtPressureCouple,
 574                 ir, md, ekind, box, x, xprime, v, f, nh_vxi, stepM);
 575         }
 576     }
 577     else
 578     {
 579         /* Use a simple and thus more efficient integration loop. */
 580         /* The simple loop does not check for particle type (so it can
 581          * be vectorized), which means we need to clear the velocities
 582          * of virtual sites in advance, when present. Note that vsite
 583          * velocities are computed after update and constraints from
 584          * their displacement.
 585          */
 586         if (md->haveVsites)
 587         {
 588             /* Note: The overhead of this loop is completely neligible */
 589             clearVsiteVelocities(start, nrend, md->ptype, v);
 590         }
 591
 592         /* We determine if we have a single T-coupling lambda value for all
 593          * atoms. That allows for better SIMD acceleration in the template.
 594          * If we do not do temperature coupling (in the run or this step),
 595          * all scaling values are 1, so we effectively have a single value.
 596          */
 597         bool haveSingleTempScaleValue = (!doTempCouple || ekind->ngtc == 1);
 598
 599         /* Extract some pointers needed by all cases */
 600         const unsigned short             *cTC           = md->cTC;
 601         gmx::ArrayRef<const t_grp_tcstat> tcstat        = ekind->tcstat;
 602         const rvec                       *invMassPerDim = md->invMassPerDim;
 603
 604         if (doParrinelloRahman)
 605         {
 606             GMX_ASSERT(!doPROffDiagonal, "updateMDLeapfrogSimple only support diagonal Parrinello-Rahman scaling matrices");
 607
 608             rvec diagM;
 609             for (int d = 0; d < DIM; d++)
 610             {
 611                 diagM[d] = M[d][d];
 612             }
 613
 614             if (haveSingleTempScaleValue)
 615             {
 616                 updateMDLeapfrogSimple
 617                 <NumTempScaleValues::single,
 618                  ApplyParrinelloRahmanVScaling::diagonal>
 619                     (start, nrend, dt, dtPressureCouple,
 620                     invMassPerDim, tcstat, cTC, diagM, x, xprime, v, f);
 621             }
 622             else
 623             {
 624                 updateMDLeapfrogSimple
 625                 <NumTempScaleValues::multiple,
 626                  ApplyParrinelloRahmanVScaling::diagonal>
 627                     (start, nrend, dt, dtPressureCouple,
 628                     invMassPerDim, tcstat, cTC, diagM, x, xprime, v, f);
 629             }
 630         }
 631         else
 632         {
 633             if (haveSingleTempScaleValue)
 634             {
 635                 /* Note that modern compilers are pretty good at vectorizing
 636                  * updateMDLeapfrogSimple(). But the SIMD version will still
 637                  * be faster because invMass lowers the cache pressure
 638                  * compared to invMassPerDim.
 639                  */
 640 #if GMX_HAVE_SIMD_UPDATE
 641                 /* Check if we can use invmass instead of invMassPerDim */
 642                 if (!md->havePartiallyFrozenAtoms)
 643                 {
 644                     updateMDLeapfrogSimpleSimd
 645                         (start, nrend, dt, md->invmass, tcstat, x, xprime, v, f);
 646                 }
 647                 else
 648 #endif
 649                 {
 650                     updateMDLeapfrogSimple
 651                     <NumTempScaleValues::single,
 652                      ApplyParrinelloRahmanVScaling::no>
 653                         (start, nrend, dt, dtPressureCouple,
 654                         invMassPerDim, tcstat, cTC, nullptr, x, xprime, v, f);
 655                 }
 656             }
 657             else
 658             {
 659                 updateMDLeapfrogSimple
 660                 <NumTempScaleValues::multiple,
 661                  ApplyParrinelloRahmanVScaling::no>
 662                     (start, nrend, dt, dtPressureCouple,
 663                     invMassPerDim, tcstat, cTC, nullptr, x, xprime, v, f);
 664             }
 665         }
 666     }
 667 }
 668
 669 static void do_update_vv_vel(int start, int nrend, real dt,
 670                              const rvec accel[], const ivec nFreeze[], const real invmass[],
 671                              const unsigned short ptype[], const unsigned short cFREEZE[],
 672                              const unsigned short cACC[], rvec v[], const rvec f[],
 673                              gmx_bool bExtended, real veta, real alpha)
 674 {
 675     int    gf = 0, ga = 0;
 676     int    n, d;
 677     real   g, mv1, mv2;
 678
 679     if (bExtended)
 680     {
 681         g        = 0.25*dt*veta*alpha;
 682         mv1      = std::exp(-g);
 683         mv2      = gmx::series_sinhx(g);
 684     }
 685     else
 686     {
 687         mv1      = 1.0;
 688         mv2      = 1.0;
 689     }
 690     for (n = start; n < nrend; n++)
 691     {
 692         real w_dt = invmass[n]*dt;
 693         if (cFREEZE)
 694         {
 695             gf   = cFREEZE[n];
 696         }
 697         if (cACC)
 698         {
 699             ga   = cACC[n];
 700         }
 701
 702         for (d = 0; d < DIM; d++)
 703         {
 704             if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
 705             {
 706                 v[n][d]             = mv1*(mv1*v[n][d] + 0.5*(w_dt*mv2*f[n][d]))+0.5*accel[ga][d]*dt;
 707             }
 708             else
 709             {
 710                 v[n][d]        = 0.0;
 711             }
 712         }
 713     }
 714 } /* do_update_vv_vel */
 715
 716 static void do_update_vv_pos(int start, int nrend, real dt,
 717                              const ivec nFreeze[],
 718                              const unsigned short ptype[], const unsigned short cFREEZE[],
 719                              const rvec x[], rvec xprime[], const rvec v[],
 720                              gmx_bool bExtended, real veta)
 721 {
 722     int    gf = 0;
 723     int    n, d;
 724     real   g, mr1, mr2;
 725
 726     /* Would it make more sense if Parrinello-Rahman was put here? */
 727     if (bExtended)
 728     {
 729         g        = 0.5*dt*veta;
 730         mr1      = std::exp(g);
 731         mr2      = gmx::series_sinhx(g);
 732     }
 733     else
 734     {
 735         mr1      = 1.0;
 736         mr2      = 1.0;
 737     }
 738
 739     for (n = start; n < nrend; n++)
 740     {
 741
 742         if (cFREEZE)
 743         {
 744             gf   = cFREEZE[n];
 745         }
 746
 747         for (d = 0; d < DIM; d++)
 748         {
 749             if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
 750             {
 751                 xprime[n][d]   = mr1*(mr1*x[n][d]+mr2*dt*v[n][d]);
 752             }
 753             else
 754             {
 755                 xprime[n][d]   = x[n][d];
 756             }
 757         }
 758     }
 759 } /* do_update_vv_pos */
 760
 761 gmx_stochd_t::gmx_stochd_t(const t_inputrec *ir)
 762 {
 763     const t_grpopts *opts = &ir->opts;
 764     const int        ngtc = opts->ngtc;
 765
 766     if (ir->eI == eiBD)
 767     {
 768         bd_rf.resize(ngtc);
 769     }
 770     else if (EI_SD(ir->eI))
 771     {
 772         sdc.resize(ngtc);
 773         sdsig.resize(ngtc);
 774
 775         for (int gt = 0; gt < ngtc; gt++)
 776         {
 777             if (opts->tau_t[gt] > 0)
 778             {
 779                 sdc[gt].em  = std::exp(-ir->delta_t/opts->tau_t[gt]);
 780             }
 781             else
 782             {
 783                 /* No friction and noise on this group */
 784                 sdc[gt].em  = 1;
 785             }
 786         }
 787     }
 788     else if (ETC_ANDERSEN(ir->etc))
 789     {
 790         randomize_group.resize(ngtc);
 791         boltzfac.resize(ngtc);
 792
 793         /* for now, assume that all groups, if randomized, are randomized at the same rate, i.e. tau_t is the same. */
 794         /* since constraint groups don't necessarily match up with temperature groups! This is checked in readir.c */
 795
 796         for (int gt = 0; gt < ngtc; gt++)
 797         {
 798             real reft = std::max<real>(0, opts->ref_t[gt]);
 799             if ((opts->tau_t[gt] > 0) && (reft > 0))  /* tau_t or ref_t = 0 means that no randomization is done */
 800             {
 801                 randomize_group[gt] = true;
 802                 boltzfac[gt]        = BOLTZ*opts->ref_t[gt];
 803             }
 804             else
 805             {
 806                 randomize_group[gt] = false;
 807             }
 808         }
 809     }
 810 }
 811
 812 void update_temperature_constants(gmx_stochd_t *sd, const t_inputrec *ir)
 813 {
 814     if (ir->eI == eiBD)
 815     {
 816         if (ir->bd_fric != 0)
 817         {
 818             for (int gt = 0; gt < ir->opts.ngtc; gt++)
 819             {
 820                 sd->bd_rf[gt] = std::sqrt(2.0*BOLTZ*ir->opts.ref_t[gt]/(ir->bd_fric*ir->delta_t));
 821             }
 822         }
 823         else
 824         {
 825             for (int gt = 0; gt < ir->opts.ngtc; gt++)
 826             {
 827                 sd->bd_rf[gt] = std::sqrt(2.0*BOLTZ*ir->opts.ref_t[gt]);
 828             }
 829         }
 830     }
 831     if (ir->eI == eiSD1)
 832     {
 833         for (int gt = 0; gt < ir->opts.ngtc; gt++)
 834         {
 835             real kT = BOLTZ*ir->opts.ref_t[gt];
 836             /* The mass is accounted for later, since this differs per atom */
 837             sd->sdsig[gt].V  = std::sqrt(kT*(1 - sd->sdc[gt].em * sd->sdc[gt].em));
 838         }
 839     }
 840 }
 841
 842 Update::Impl::Impl(const t_inputrec    *ir, BoxDeformation *boxDeformation)
 843 {
 844     sd = std::make_unique<gmx_stochd_t>(ir);
 845     update_temperature_constants(sd.get(), ir);
 846     xp.resizeWithPadding(0);
 847     deform = boxDeformation;
 848 }
 849
 850 void Update::setNumAtoms(int nAtoms)
 851 {
 852
 853     impl_->xp.resizeWithPadding(nAtoms);
 854 }
 855
 856 /*! \brief Sets the SD update type */
 857 enum class SDUpdate : int
 858 {
 859     ForcesOnly, FrictionAndNoiseOnly, Combined
 860 };
 861
 862 /*! \brief SD integrator update
 863  *
 864  * Two phases are required in the general case of a constrained
 865  * update, the first phase from the contribution of forces, before
 866  * applying constraints, and then a second phase applying the friction
 867  * and noise, and then further constraining. For details, see
 868  * Goga2012.
 869  *
 870  * Without constraints, the two phases can be combined, for
 871  * efficiency.
 872  *
 873  * Thus three instantiations of this templated function will be made,
 874  * two with only one contribution, and one with both contributions. */
 875 template <SDUpdate updateType>
 876 static void
 877 doSDUpdateGeneral(const gmx_stochd_t &sd,
 878                   int start, int nrend, real dt,
 879                   const rvec accel[], const ivec nFreeze[],
 880                   const real invmass[], const unsigned short ptype[],
 881                   const unsigned short cFREEZE[], const unsigned short cACC[],
 882                   const unsigned short cTC[],
 883                   const rvec x[], rvec xprime[], rvec v[], const rvec f[],
 884                   int64_t step, int seed, const int *gatindex)
 885 {
 886     // cTC, cACC and cFREEZE can be nullptr any time, but various
 887     // instantiations do not make sense with particular pointer
 888     // values.
 889     if (updateType == SDUpdate::ForcesOnly)
 890     {
 891         GMX_ASSERT(f != nullptr, "SD update with only forces requires forces");
 892         GMX_ASSERT(cTC == nullptr, "SD update with only forces cannot handle temperature groups");
 893     }
 894     if (updateType == SDUpdate::FrictionAndNoiseOnly)
 895     {
 896         GMX_ASSERT(f == nullptr, "SD update with only noise cannot handle forces");
 897         GMX_ASSERT(cACC == nullptr, "SD update with only noise cannot handle acceleration groups");
 898     }
 899     if (updateType == SDUpdate::Combined)
 900     {
 901         GMX_ASSERT(f != nullptr, "SD update with forces and noise requires forces");
 902     }
 903
 904     // Even 0 bits internal counter gives 2x64 ints (more than enough for three table lookups)
 905     gmx::ThreeFry2x64<0> rng(seed, gmx::RandomDomain::UpdateCoordinates);
 906     gmx::TabulatedNormalDistribution<real, 14> dist;
 907
 908     for (int n = start; n < nrend; n++)
 909     {
 910         int globalAtomIndex = gatindex ? gatindex[n] : n;
 911         rng.restart(step, globalAtomIndex);
 912         dist.reset();
 913
 914         real inverseMass     = invmass[n];
 915         real invsqrtMass     = std::sqrt(inverseMass);
 916
 917         int  freezeGroup       = cFREEZE ? cFREEZE[n] : 0;
 918         int  accelerationGroup = cACC ? cACC[n] : 0;
 919         int  temperatureGroup  = cTC ? cTC[n] : 0;
 920
 921         for (int d = 0; d < DIM; d++)
 922         {
 923             if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[freezeGroup][d])
 924             {
 925                 if (updateType == SDUpdate::ForcesOnly)
 926                 {
 927                     real vn      = v[n][d] + (inverseMass*f[n][d] + accel[accelerationGroup][d])*dt;
 928                     v[n][d]      = vn;
 929                     // Simple position update.
 930                     xprime[n][d] = x[n][d] + v[n][d]*dt;
 931                 }
 932                 else if (updateType == SDUpdate::FrictionAndNoiseOnly)
 933                 {
 934                     real vn      = v[n][d];
 935                     v[n][d]      = (vn*sd.sdc[temperatureGroup].em +
 936                                     invsqrtMass*sd.sdsig[temperatureGroup].V*dist(rng));
 937                     // The previous phase already updated the
 938                     // positions with a full v*dt term that must
 939                     // now be half removed.
 940                     xprime[n][d] = xprime[n][d] + 0.5*(v[n][d] - vn)*dt;
 941                 }
 942                 else
 943                 {
 944                     real vn      = v[n][d] + (inverseMass*f[n][d] + accel[accelerationGroup][d])*dt;
 945                     v[n][d]      = (vn*sd.sdc[temperatureGroup].em +
 946                                     invsqrtMass*sd.sdsig[temperatureGroup].V*dist(rng));
 947                     // Here we include half of the friction+noise
 948                     // update of v into the position update.
 949                     xprime[n][d] = x[n][d] + 0.5*(vn + v[n][d])*dt;
 950                 }
 951             }
 952             else
 953             {
 954                 // When using constraints, the update is split into
 955                 // two phases, but we only need to zero the update of
 956                 // virtual, shell or frozen particles in at most one
 957                 // of the phases.
 958                 if (updateType != SDUpdate::FrictionAndNoiseOnly)
 959                 {
 960                     v[n][d]      = 0.0;
 961                     xprime[n][d] = x[n][d];
 962                 }
 963             }
 964         }
 965     }
 966 }
 967
 968 static void do_update_bd(int start, int nrend, real dt,
 969                          const ivec nFreeze[],
 970                          const real invmass[], const unsigned short ptype[],
 971                          const unsigned short cFREEZE[], const unsigned short cTC[],
 972                          const rvec x[], rvec xprime[], rvec v[],
 973                          const rvec f[], real friction_coefficient,
 974                          const real *rf, int64_t step, int seed,
 975                          const int* gatindex)
 976 {
 977     /* note -- these appear to be full step velocities . . .  */
 978     int    gf = 0, gt = 0;
 979     real   vn;
 980     real   invfr = 0;
 981     int    n, d;
 982     // Use 1 bit of internal counters to give us 2*2 64-bits values per stream
 983     // Each 64-bit value is enough for 4 normal distribution table numbers.
 984     gmx::ThreeFry2x64<0> rng(seed, gmx::RandomDomain::UpdateCoordinates);
 985     gmx::TabulatedNormalDistribution<real, 14> dist;
 986
 987     if (friction_coefficient != 0)
 988     {
 989         invfr = 1.0/friction_coefficient;
 990     }
 991
 992     for (n = start; (n < nrend); n++)
 993     {
 994         int  ng  = gatindex ? gatindex[n] : n;
 995
 996         rng.restart(step, ng);
 997         dist.reset();
 998
 999         if (cFREEZE)
1000         {
1001             gf = cFREEZE[n];
1002         }
1003         if (cTC)
1004         {
1005             gt = cTC[n];
1006         }
1007         for (d = 0; (d < DIM); d++)
1008         {
1009             if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
1010             {
1011                 if (friction_coefficient != 0)
1012                 {
1013                     vn = invfr*f[n][d] + rf[gt]*dist(rng);
1014                 }
1015                 else
1016                 {
1017                     /* NOTE: invmass = 2/(mass*friction_constant*dt) */
1018                     vn = 0.5*invmass[n]*f[n][d]*dt
1019                         + std::sqrt(0.5*invmass[n])*rf[gt]*dist(rng);
1020                 }
1021
1022                 v[n][d]      = vn;
1023                 xprime[n][d] = x[n][d]+vn*dt;
1024             }
1025             else
1026             {
1027                 v[n][d]      = 0.0;
1028                 xprime[n][d] = x[n][d];
1029             }
1030         }
1031     }
1032 }
1033
1034 static void calc_ke_part_normal(const rvec v[], const t_grpopts *opts, const t_mdatoms *md,
1035                                 gmx_ekindata_t *ekind, t_nrnb *nrnb, gmx_bool bEkinAveVel)
1036 {
1037     int                         g;
1038     gmx::ArrayRef<t_grp_tcstat> tcstat  = ekind->tcstat;
1039     gmx::ArrayRef<t_grp_acc>    grpstat = ekind->grpstat;
1040
1041     /* three main: VV with AveVel, vv with AveEkin, leap with AveEkin.  Leap with AveVel is also
1042        an option, but not supported now.
1043        bEkinAveVel: If TRUE, we sum into ekin, if FALSE, into ekinh.
1044      */
1045
1046     /* group velocities are calculated in update_ekindata and
1047      * accumulated in acumulate_groups.
1048      * Now the partial global and groups ekin.
1049      */
1050     for (g = 0; (g < opts->ngtc); g++)
1051     {
1052         copy_mat(tcstat[g].ekinh, tcstat[g].ekinh_old);
1053         if (bEkinAveVel)
1054         {
1055             clear_mat(tcstat[g].ekinf);
1056             tcstat[g].ekinscalef_nhc = 1.0;   /* need to clear this -- logic is complicated! */
1057         }
1058         else
1059         {
1060             clear_mat(tcstat[g].ekinh);
1061         }
1062     }
1063     ekind->dekindl_old = ekind->dekindl;
1064     int nthread        = gmx_omp_nthreads_get(emntUpdate);
1065
1066 #pragma omp parallel for num_threads(nthread) schedule(static)
1067     for (int thread = 0; thread < nthread; thread++)
1068     {
1069         // This OpenMP only loops over arrays and does not call any functions
1070         // or memory allocation. It should not be able to throw, so for now
1071         // we do not need a try/catch wrapper.
1072         int     start_t, end_t, n;
1073         int     ga, gt;
1074         rvec    v_corrt;
1075         real    hm;
1076         int     d, m;
1077         matrix *ekin_sum;
1078         real   *dekindl_sum;
1079
1080         start_t = ((thread+0)*md->homenr)/nthread;
1081         end_t   = ((thread+1)*md->homenr)/nthread;
1082
1083         ekin_sum    = ekind->ekin_work[thread];
1084         dekindl_sum = ekind->dekindl_work[thread];
1085
1086         for (gt = 0; gt < opts->ngtc; gt++)
1087         {
1088             clear_mat(ekin_sum[gt]);
1089         }
1090         *dekindl_sum = 0.0;
1091
1092         ga = 0;
1093         gt = 0;
1094         for (n = start_t; n < end_t; n++)
1095         {
1096             if (md->cACC)
1097             {
1098                 ga = md->cACC[n];
1099             }
1100             if (md->cTC)
1101             {
1102                 gt = md->cTC[n];
1103             }
1104             hm   = 0.5*md->massT[n];
1105
1106             for (d = 0; (d < DIM); d++)
1107             {
1108                 v_corrt[d]  = v[n][d]  - grpstat[ga].u[d];
1109             }
1110             for (d = 0; (d < DIM); d++)
1111             {
1112                 for (m = 0; (m < DIM); m++)
1113                 {
1114                     /* if we're computing a full step velocity, v_corrt[d] has v(t).  Otherwise, v(t+dt/2) */
1115                     ekin_sum[gt][m][d] += hm*v_corrt[m]*v_corrt[d];
1116                 }
1117             }
1118             if (md->nMassPerturbed && md->bPerturbed[n])
1119             {
1120                 *dekindl_sum +=
1121                     0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
1122             }
1123         }
1124     }
1125
1126     ekind->dekindl = 0;
1127     for (int thread = 0; thread < nthread; thread++)
1128     {
1129         for (g = 0; g < opts->ngtc; g++)
1130         {
1131             if (bEkinAveVel)
1132             {
1133                 m_add(tcstat[g].ekinf, ekind->ekin_work[thread][g],
1134                       tcstat[g].ekinf);
1135             }
1136             else
1137             {
1138                 m_add(tcstat[g].ekinh, ekind->ekin_work[thread][g],
1139                       tcstat[g].ekinh);
1140             }
1141         }
1142
1143         ekind->dekindl += *ekind->dekindl_work[thread];
1144     }
1145
1146     inc_nrnb(nrnb, eNR_EKIN, md->homenr);
1147 }
1148
1149 static void calc_ke_part_visc(const matrix box, const rvec x[], const rvec v[],
1150                               const t_grpopts *opts, const t_mdatoms *md,
1151                               gmx_ekindata_t *ekind,
1152                               t_nrnb *nrnb, gmx_bool bEkinAveVel)
1153 {
1154     int                         start = 0, homenr = md->homenr;
1155     int                         g, d, n, m, gt = 0;
1156     rvec                        v_corrt;
1157     real                        hm;
1158     gmx::ArrayRef<t_grp_tcstat> tcstat = ekind->tcstat;
1159     t_cos_acc                  *cosacc = &(ekind->cosacc);
1160     real                        dekindl;
1161     real                        fac, cosz;
1162     double                      mvcos;
1163
1164     for (g = 0; g < opts->ngtc; g++)
1165     {
1166         copy_mat(ekind->tcstat[g].ekinh, ekind->tcstat[g].ekinh_old);
1167         clear_mat(ekind->tcstat[g].ekinh);
1168     }
1169     ekind->dekindl_old = ekind->dekindl;
1170
1171     fac     = 2*M_PI/box[ZZ][ZZ];
1172     mvcos   = 0;
1173     dekindl = 0;
1174     for (n = start; n < start+homenr; n++)
1175     {
1176         if (md->cTC)
1177         {
1178             gt = md->cTC[n];
1179         }
1180         hm   = 0.5*md->massT[n];
1181
1182         /* Note that the times of x and v differ by half a step */
1183         /* MRS -- would have to be changed for VV */
1184         cosz         = std::cos(fac*x[n][ZZ]);
1185         /* Calculate the amplitude of the new velocity profile */
1186         mvcos       += 2*cosz*md->massT[n]*v[n][XX];
1187
1188         copy_rvec(v[n], v_corrt);
1189         /* Subtract the profile for the kinetic energy */
1190         v_corrt[XX] -= cosz*cosacc->vcos;
1191         for (d = 0; (d < DIM); d++)
1192         {
1193             for (m = 0; (m < DIM); m++)
1194             {
1195                 /* if we're computing a full step velocity, v_corrt[d] has v(t).  Otherwise, v(t+dt/2) */
1196                 if (bEkinAveVel)
1197                 {
1198                     tcstat[gt].ekinf[m][d] += hm*v_corrt[m]*v_corrt[d];
1199                 }
1200                 else
1201                 {
1202                     tcstat[gt].ekinh[m][d] += hm*v_corrt[m]*v_corrt[d];
1203                 }
1204             }
1205         }
1206         if (md->nPerturbed && md->bPerturbed[n])
1207         {
1208             /* The minus sign here might be confusing.
1209              * The kinetic contribution from dH/dl doesn't come from
1210              * d m(l)/2 v^2 / dl, but rather from d p^2/2m(l) / dl,
1211              * where p are the momenta. The difference is only a minus sign.
1212              */
1213             dekindl -= 0.5*(md->massB[n] - md->massA[n])*iprod(v_corrt, v_corrt);
1214         }
1215     }
1216     ekind->dekindl = dekindl;
1217     cosacc->mvcos  = mvcos;
1218
1219     inc_nrnb(nrnb, eNR_EKIN, homenr);
1220 }
1221
1222 void calc_ke_part(
1223         const rvec *x, const rvec *v, const matrix box,
1224         const t_grpopts *opts, const t_mdatoms *md,
1225         gmx_ekindata_t *ekind, t_nrnb *nrnb, gmx_bool bEkinAveVel)
1226 {
1227     if (ekind->cosacc.cos_accel == 0)
1228     {
1229         calc_ke_part_normal(v, opts, md, ekind, nrnb, bEkinAveVel);
1230     }
1231     else
1232     {
1233         calc_ke_part_visc(box, x, v, opts, md, ekind, nrnb, bEkinAveVel);
1234     }
1235 }
1236
1237 extern void init_ekinstate(ekinstate_t *ekinstate, const t_inputrec *ir)
1238 {
1239     ekinstate->ekin_n = ir->opts.ngtc;
1240     snew(ekinstate->ekinh, ekinstate->ekin_n);
1241     snew(ekinstate->ekinf, ekinstate->ekin_n);
1242     snew(ekinstate->ekinh_old, ekinstate->ekin_n);
1243     ekinstate->ekinscalef_nhc.resize(ekinstate->ekin_n);
1244     ekinstate->ekinscaleh_nhc.resize(ekinstate->ekin_n);
1245     ekinstate->vscale_nhc.resize(ekinstate->ekin_n);
1246     ekinstate->dekindl          = 0;
1247     ekinstate->mvcos            = 0;
1248     ekinstate->hasReadEkinState = false;
1249 }
1250
1251 void update_ekinstate(ekinstate_t *ekinstate, const gmx_ekindata_t *ekind)
1252 {
1253     int i;
1254
1255     for (i = 0; i < ekinstate->ekin_n; i++)
1256     {
1257         copy_mat(ekind->tcstat[i].ekinh, ekinstate->ekinh[i]);
1258         copy_mat(ekind->tcstat[i].ekinf, ekinstate->ekinf[i]);
1259         copy_mat(ekind->tcstat[i].ekinh_old, ekinstate->ekinh_old[i]);
1260         ekinstate->ekinscalef_nhc[i] = ekind->tcstat[i].ekinscalef_nhc;
1261         ekinstate->ekinscaleh_nhc[i] = ekind->tcstat[i].ekinscaleh_nhc;
1262         ekinstate->vscale_nhc[i]     = ekind->tcstat[i].vscale_nhc;
1263     }
1264
1265     copy_mat(ekind->ekin, ekinstate->ekin_total);
1266     ekinstate->dekindl = ekind->dekindl;
1267     ekinstate->mvcos   = ekind->cosacc.mvcos;
1268
1269 }
1270
1271 void restore_ekinstate_from_state(const t_commrec *cr,
1272                                   gmx_ekindata_t *ekind, const ekinstate_t *ekinstate)
1273 {
1274     int i, n;
1275
1276     if (MASTER(cr))
1277     {
1278         for (i = 0; i < ekinstate->ekin_n; i++)
1279         {
1280             copy_mat(ekinstate->ekinh[i], ekind->tcstat[i].ekinh);
1281             copy_mat(ekinstate->ekinf[i], ekind->tcstat[i].ekinf);
1282             copy_mat(ekinstate->ekinh_old[i], ekind->tcstat[i].ekinh_old);
1283             ekind->tcstat[i].ekinscalef_nhc = ekinstate->ekinscalef_nhc[i];
1284             ekind->tcstat[i].ekinscaleh_nhc = ekinstate->ekinscaleh_nhc[i];
1285             ekind->tcstat[i].vscale_nhc     = ekinstate->vscale_nhc[i];
1286         }
1287
1288         copy_mat(ekinstate->ekin_total, ekind->ekin);
1289
1290         ekind->dekindl      = ekinstate->dekindl;
1291         ekind->cosacc.mvcos = ekinstate->mvcos;
1292         n                   = ekinstate->ekin_n;
1293     }
1294
1295     if (PAR(cr))
1296     {
1297         gmx_bcast(sizeof(n), &n, cr);
1298         for (i = 0; i < n; i++)
1299         {
1300             gmx_bcast(DIM*DIM*sizeof(ekind->tcstat[i].ekinh[0][0]),
1301                       ekind->tcstat[i].ekinh[0], cr);
1302             gmx_bcast(DIM*DIM*sizeof(ekind->tcstat[i].ekinf[0][0]),
1303                       ekind->tcstat[i].ekinf[0], cr);
1304             gmx_bcast(DIM*DIM*sizeof(ekind->tcstat[i].ekinh_old[0][0]),
1305                       ekind->tcstat[i].ekinh_old[0], cr);
1306
1307             gmx_bcast(sizeof(ekind->tcstat[i].ekinscalef_nhc),
1308                       &(ekind->tcstat[i].ekinscalef_nhc), cr);
1309             gmx_bcast(sizeof(ekind->tcstat[i].ekinscaleh_nhc),
1310                       &(ekind->tcstat[i].ekinscaleh_nhc), cr);
1311             gmx_bcast(sizeof(ekind->tcstat[i].vscale_nhc),
1312                       &(ekind->tcstat[i].vscale_nhc), cr);
1313         }
1314         gmx_bcast(DIM*DIM*sizeof(ekind->ekin[0][0]),
1315                   ekind->ekin[0], cr);
1316
1317         gmx_bcast(sizeof(ekind->dekindl), &ekind->dekindl, cr);
1318         gmx_bcast(sizeof(ekind->cosacc.mvcos), &ekind->cosacc.mvcos, cr);
1319     }
1320 }
1321
1322 void update_tcouple(int64_t           step,
1323                     const t_inputrec *inputrec,
1324                     t_state          *state,
1325                     gmx_ekindata_t   *ekind,
1326                     const t_extmass  *MassQ,
1327                     const t_mdatoms  *md)
1328
1329 {
1330     // This condition was explicitly checked in previous version, but should have never been satisfied
1331     GMX_ASSERT(!(EI_VV(inputrec->eI) &&
1332                  (inputrecNvtTrotter(inputrec) || inputrecNptTrotter(inputrec) || inputrecNphTrotter(inputrec))),
1333                "Temperature coupling was requested with velocity verlet and trotter");
1334
1335     bool doTemperatureCoupling = false;
1336
1337     // For VV temperature coupling parameters are updated on the current
1338     // step, for the others - one step before.
1339     if (inputrec->etc == etcNO)
1340     {
1341         doTemperatureCoupling = false;
1342     }
1343     else if (EI_VV(inputrec->eI))
1344     {
1345         doTemperatureCoupling = do_per_step(step, inputrec->nsttcouple);
1346     }
1347     else
1348     {
1349         doTemperatureCoupling = do_per_step(step + inputrec->nsttcouple - 1, inputrec->nsttcouple);
1350     }
1351
1352     if (doTemperatureCoupling)
1353     {
1354         real dttc = inputrec->nsttcouple*inputrec->delta_t;
1355
1356         //TODO: berendsen_tcoupl(...), nosehoover_tcoupl(...) and vrescale_tcoupl(...) update
1357         //      temperature coupling parameters, which should be reflected in the name of these
1358         //      subroutines
1359         switch (inputrec->etc)
1360         {
1361             case etcNO:
1362                 break;
1363             case etcBERENDSEN:
1364                 berendsen_tcoupl(inputrec, ekind, dttc, state->therm_integral);
1365                 break;
1366             case etcNOSEHOOVER:
1367                 nosehoover_tcoupl(&(inputrec->opts), ekind, dttc,
1368                                   state->nosehoover_xi.data(), state->nosehoover_vxi.data(), MassQ);
1369                 break;
1370             case etcVRESCALE:
1371                 vrescale_tcoupl(inputrec, step, ekind, dttc,
1372                                 state->therm_integral.data());
1373                 break;
1374         }
1375         /* rescale in place here */
1376         if (EI_VV(inputrec->eI))
1377         {
1378             rescale_velocities(ekind, md, 0, md->homenr, state->v.rvec_array());
1379         }
1380     }
1381     else
1382     {
1383         // Set the T scaling lambda to 1 to have no scaling
1384         // TODO: Do we have to do it on every non-t-couple step?
1385         for (int i = 0; (i < inputrec->opts.ngtc); i++)
1386         {
1387             ekind->tcstat[i].lambda = 1.0;
1388         }
1389     }
1390 }
1391
1392 void getThreadAtomRange(int numThreads, int threadIndex, int numAtoms,
1393                         int *startAtom, int *endAtom)
1394 {
1395 #if GMX_HAVE_SIMD_UPDATE
1396     constexpr int blockSize = GMX_SIMD_REAL_WIDTH;
1397 #else
1398     constexpr int blockSize = 1;
1399 #endif
1400     int           numBlocks = (numAtoms + blockSize - 1)/blockSize;
1401
1402     *startAtom    = ((numBlocks* threadIndex     )/numThreads)*blockSize;
1403     *endAtom      = ((numBlocks*(threadIndex + 1))/numThreads)*blockSize;
1404     if (threadIndex == numThreads - 1)
1405     {
1406         *endAtom  = numAtoms;
1407     }
1408 }
1409
1410 void update_pcouple_before_coordinates(FILE             *fplog,
1411                                        int64_t           step,
1412                                        const t_inputrec *inputrec,
1413                                        t_state          *state,
1414                                        matrix            parrinellorahmanMu,
1415                                        matrix            M,
1416                                        gmx_bool          bInitStep)
1417 {
1418     /* Berendsen P-coupling is completely handled after the coordinate update.
1419      * Trotter P-coupling is handled by separate calls to trotter_update().
1420      */
1421     if (inputrec->epc == epcPARRINELLORAHMAN &&
1422         do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
1423     {
1424         real dtpc = inputrec->nstpcouple*inputrec->delta_t;
1425
1426         parrinellorahman_pcoupl(fplog, step, inputrec, dtpc, state->pres_prev,
1427                                 state->box, state->box_rel, state->boxv,
1428                                 M, parrinellorahmanMu, bInitStep);
1429     }
1430 }
1431
1432 void constrain_velocities(int64_t                        step,
1433                           real                          *dvdlambda, /* the contribution to be added to the bonded interactions */
1434                           t_state                       *state,
1435                           tensor                         vir_part,
1436                           gmx::Constraints              *constr,
1437                           gmx_bool                       bCalcVir,
1438                           bool                           do_log,
1439                           bool                           do_ene)
1440 {
1441     if (!constr)
1442     {
1443         return;
1444     }
1445
1446     /*
1447      *  Steps (7C, 8C)
1448      *  APPLY CONSTRAINTS:
1449      *  BLOCK SHAKE
1450      */
1451
1452     {
1453         tensor vir_con;
1454
1455         /* clear out constraints before applying */
1456         clear_mat(vir_part);
1457
1458         /* Constrain the coordinates upd->xp */
1459         constr->apply(do_log, do_ene,
1460                       step, 1, 1.0,
1461                       state->x.rvec_array(), state->v.rvec_array(), state->v.rvec_array(),
1462                       state->box,
1463                       state->lambda[efptBONDED], dvdlambda,
1464                       nullptr, bCalcVir ? &vir_con : nullptr, ConstraintVariable::Velocities);
1465
1466         if (bCalcVir)
1467         {
1468             m_add(vir_part, vir_con, vir_part);
1469         }
1470     }
1471 }
1472
1473 void constrain_coordinates(int64_t           step,
1474                            real             *dvdlambda, /* the contribution to be added to the bonded interactions */
1475                            t_state          *state,
1476                            tensor            vir_part,
1477                            Update           *upd,
1478                            gmx::Constraints *constr,
1479                            gmx_bool          bCalcVir,
1480                            bool              do_log,
1481                            bool              do_ene)
1482 {
1483     if (!constr)
1484     {
1485         return;
1486     }
1487
1488     {
1489         tensor vir_con;
1490
1491         /* clear out constraints before applying */
1492         clear_mat(vir_part);
1493
1494         /* Constrain the coordinates upd->xp */
1495         constr->apply(do_log, do_ene,
1496                       step, 1, 1.0,
1497                       state->x.rvec_array(), upd->xp()->rvec_array(), nullptr,
1498                       state->box,
1499                       state->lambda[efptBONDED], dvdlambda,
1500                       as_rvec_array(state->v.data()), bCalcVir ? &vir_con : nullptr, ConstraintVariable::Positions);
1501
1502         if (bCalcVir)
1503         {
1504             m_add(vir_part, vir_con, vir_part);
1505         }
1506     }
1507 }
1508
1509 void
1510 update_sd_second_half(int64_t           step,
1511                       real             *dvdlambda, /* the contribution to be added to the bonded interactions */
1512                       const t_inputrec *inputrec,  /* input record and box stuff        */
1513                       const t_mdatoms  *md,
1514                       t_state          *state,
1515                       const t_commrec  *cr,
1516                       t_nrnb           *nrnb,
1517                       gmx_wallcycle_t   wcycle,
1518                       Update           *upd,
1519                       gmx::Constraints *constr,
1520                       bool              do_log,
1521                       bool              do_ene)
1522 {
1523     if (!constr)
1524     {
1525         return;
1526     }
1527     if (inputrec->eI == eiSD1)
1528     {
1529         int homenr = md->homenr;
1530
1531         /* Cast delta_t from double to real to make the integrators faster.
1532          * The only reason for having delta_t double is to get accurate values
1533          * for t=delta_t*step when step is larger than float precision.
1534          * For integration dt the accuracy of real suffices, since with
1535          * integral += dt*integrand the increment is nearly always (much) smaller
1536          * than the integral (and the integrand has real precision).
1537          */
1538         real dt     = inputrec->delta_t;
1539
1540         wallcycle_start(wcycle, ewcUPDATE);
1541
1542         int nth = gmx_omp_nthreads_get(emntUpdate);
1543
1544 #pragma omp parallel for num_threads(nth) schedule(static)
1545         for (int th = 0; th < nth; th++)
1546         {
1547             try
1548             {
1549                 int start_th, end_th;
1550                 getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
1551
1552                 doSDUpdateGeneral<SDUpdate::FrictionAndNoiseOnly>
1553                     (*upd->sd(),
1554                     start_th, end_th, dt,
1555                     inputrec->opts.acc, inputrec->opts.nFreeze,
1556                     md->invmass, md->ptype,
1557                     md->cFREEZE, nullptr, md->cTC,
1558                     state->x.rvec_array(), upd->xp()->rvec_array(),
1559                     state->v.rvec_array(), nullptr,
1560                     step, inputrec->ld_seed,
1561                     DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1562             }
1563             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1564         }
1565         inc_nrnb(nrnb, eNR_UPDATE, homenr);
1566         wallcycle_stop(wcycle, ewcUPDATE);
1567
1568         /* Constrain the coordinates upd->xp for half a time step */
1569         constr->apply(do_log, do_ene,
1570                       step, 1, 0.5,
1571                       state->x.rvec_array(), upd->xp()->rvec_array(), nullptr,
1572                       state->box,
1573                       state->lambda[efptBONDED], dvdlambda,
1574                       as_rvec_array(state->v.data()), nullptr, ConstraintVariable::Positions);
1575     }
1576 }
1577
1578 void finish_update(const t_inputrec       *inputrec, /* input record and box stuff      */
1579                    const t_mdatoms        *md,
1580                    t_state                *state,
1581                    const t_graph          *graph,
1582                    t_nrnb                 *nrnb,
1583                    gmx_wallcycle_t         wcycle,
1584                    Update                 *upd,
1585                    const gmx::Constraints *constr)
1586 {
1587     int homenr = md->homenr;
1588
1589     /* We must always unshift after updating coordinates; if we did not shake
1590        x was shifted in do_force */
1591
1592     /* NOTE Currently we always integrate to a temporary buffer and
1593      * then copy the results back. */
1594     {
1595         wallcycle_start_nocount(wcycle, ewcUPDATE);
1596
1597         if (md->cFREEZE != nullptr && constr != nullptr)
1598         {
1599             /* If we have atoms that are frozen along some, but not all
1600              * dimensions, then any constraints will have moved them also along
1601              * the frozen dimensions. To freeze such degrees of freedom
1602              * we copy them back here to later copy them forward. It would
1603              * be more elegant and slightly more efficient to copies zero
1604              * times instead of twice, but the graph case below prevents this.
1605              */
1606             const ivec *nFreeze                     = inputrec->opts.nFreeze;
1607             bool        partialFreezeAndConstraints = false;
1608             for (int g = 0; g < inputrec->opts.ngfrz; g++)
1609             {
1610                 int numFreezeDim = nFreeze[g][XX] + nFreeze[g][YY] + nFreeze[g][ZZ];
1611                 if (numFreezeDim > 0 && numFreezeDim < 3)
1612                 {
1613                     partialFreezeAndConstraints = true;
1614                 }
1615             }
1616             if (partialFreezeAndConstraints)
1617             {
1618                 auto xp = makeArrayRef(*upd->xp()).subArray(0, homenr);
1619                 auto x  = makeConstArrayRef(state->x).subArray(0, homenr);
1620                 for (int i = 0; i < homenr; i++)
1621                 {
1622                     int g = md->cFREEZE[i];
1623
1624                     for (int d = 0; d < DIM; d++)
1625                     {
1626                         if (nFreeze[g][d])
1627                         {
1628                             xp[i][d] = x[i][d];
1629                         }
1630                     }
1631                 }
1632             }
1633         }
1634
1635         if (graph && (graph->nnodes > 0))
1636         {
1637             unshift_x(graph, state->box, state->x.rvec_array(), upd->xp()->rvec_array());
1638             if (TRICLINIC(state->box))
1639             {
1640                 inc_nrnb(nrnb, eNR_SHIFTX, 2*graph->nnodes);
1641             }
1642             else
1643             {
1644                 inc_nrnb(nrnb, eNR_SHIFTX, graph->nnodes);
1645             }
1646         }
1647         else
1648         {
1649             auto xp  = makeConstArrayRef(*upd->xp()).subArray(0, homenr);
1650             auto x   = makeArrayRef(state->x).subArray(0, homenr);
1651
1652
1653             int gmx_unused nth = gmx_omp_nthreads_get(emntUpdate);
1654 #pragma omp parallel for num_threads(nth) schedule(static)
1655             for (int i = 0; i < homenr; i++)
1656             {
1657                 // Trivial statement, does not throw
1658                 x[i] = xp[i];
1659             }
1660         }
1661         wallcycle_stop(wcycle, ewcUPDATE);
1662     }
1663     /* ############# END the update of velocities and positions ######### */
1664 }
1665
1666 void update_pcouple_after_coordinates(FILE             *fplog,
1667                                       int64_t           step,
1668                                       const t_inputrec *inputrec,
1669                                       const t_mdatoms  *md,
1670                                       const matrix      pressure,
1671                                       const matrix      forceVirial,
1672                                       const matrix      constraintVirial,
1673                                       const matrix      parrinellorahmanMu,
1674                                       t_state          *state,
1675                                       t_nrnb           *nrnb,
1676                                       Update           *upd)
1677 {
1678     int  start  = 0;
1679     int  homenr = md->homenr;
1680
1681     /* Cast to real for faster code, no loss in precision (see comment above) */
1682     real dt     = inputrec->delta_t;
1683
1684
1685     /* now update boxes */
1686     switch (inputrec->epc)
1687     {
1688         case (epcNO):
1689             break;
1690         case (epcBERENDSEN):
1691             if (do_per_step(step, inputrec->nstpcouple))
1692             {
1693                 real   dtpc = inputrec->nstpcouple*dt;
1694                 matrix mu;
1695                 berendsen_pcoupl(fplog, step, inputrec, dtpc,
1696                                  pressure, state->box,
1697                                  forceVirial, constraintVirial,
1698                                  mu, &state->baros_integral);
1699                 berendsen_pscale(inputrec, mu, state->box, state->box_rel,
1700                                  start, homenr, state->x.rvec_array(),
1701                                  md->cFREEZE, nrnb);
1702             }
1703             break;
1704         case (epcPARRINELLORAHMAN):
1705             if (do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
1706             {
1707                 /* The box velocities were updated in do_pr_pcoupl,
1708                  * but we dont change the box vectors until we get here
1709                  * since we need to be able to shift/unshift above.
1710                  */
1711                 real dtpc = inputrec->nstpcouple*dt;
1712                 for (int i = 0; i < DIM; i++)
1713                 {
1714                     for (int m = 0; m <= i; m++)
1715                     {
1716                         state->box[i][m] += dtpc*state->boxv[i][m];
1717                     }
1718                 }
1719                 preserve_box_shape(inputrec, state->box_rel, state->box);
1720
1721                 /* Scale the coordinates */
1722                 auto x = state->x.rvec_array();
1723                 for (int n = start; n < start + homenr; n++)
1724                 {
1725                     tmvmul_ur0(parrinellorahmanMu, x[n], x[n]);
1726                 }
1727             }
1728             break;
1729         case (epcMTTK):
1730             switch (inputrec->epct)
1731             {
1732                 case (epctISOTROPIC):
1733                     /* DIM * eta = ln V.  so DIM*eta_new = DIM*eta_old + DIM*dt*veta =>
1734                        ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) =>
1735                        Side length scales as exp(veta*dt) */
1736
1737                     msmul(state->box, std::exp(state->veta*dt), state->box);
1738
1739                     /* Relate veta to boxv.  veta = d(eta)/dT = (1/DIM)*1/V dV/dT.
1740                        o               If we assume isotropic scaling, and box length scaling
1741                        factor L, then V = L^DIM (det(M)).  So dV/dt = DIM
1742                        L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt.  The
1743                        determinant of B is L^DIM det(M), and the determinant
1744                        of dB/dt is (dL/dT)^DIM det (M).  veta will be
1745                        (det(dB/dT)/det(B))^(1/3).  Then since M =
1746                        B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */
1747
1748                     msmul(state->box, state->veta, state->boxv);
1749                     break;
1750                 default:
1751                     break;
1752             }
1753             break;
1754         default:
1755             break;
1756     }
1757
1758     if (upd->deform())
1759     {
1760         auto localX = makeArrayRef(state->x).subArray(start, homenr);
1761         upd->deform()->apply(localX, state->box, step);
1762     }
1763 }
1764
1765 void update_coords(int64_t                                   step,
1766                    const t_inputrec                         *inputrec, /* input record and box stuff    */
1767                    const t_mdatoms                          *md,
1768                    t_state                                  *state,
1769                    gmx::ArrayRefWithPadding<const gmx::RVec> f,
1770                    const t_fcdata                           *fcd,
1771                    const gmx_ekindata_t                     *ekind,
1772                    const matrix                              M,
1773                    Update                                   *upd,
1774                    int                                       UpdatePart,
1775                    const t_commrec                          *cr, /* these shouldn't be here -- need to think about it */
1776                    const gmx::Constraints                   *constr)
1777 {
1778     gmx_bool bDoConstr = (nullptr != constr);
1779
1780     /* Running the velocity half does nothing except for velocity verlet */
1781     if ((UpdatePart == etrtVELOCITY1 || UpdatePart == etrtVELOCITY2) &&
1782         !EI_VV(inputrec->eI))
1783     {
1784         gmx_incons("update_coords called for velocity without VV integrator");
1785     }
1786
1787     int  homenr = md->homenr;
1788
1789     /* Cast to real for faster code, no loss in precision (see comment above) */
1790     real dt     = inputrec->delta_t;
1791
1792     /* We need to update the NMR restraint history when time averaging is used */
1793     if (state->flags & (1<<estDISRE_RM3TAV))
1794     {
1795         update_disres_history(fcd, &state->hist);
1796     }
1797     if (state->flags & (1<<estORIRE_DTAV))
1798     {
1799         update_orires_history(fcd, &state->hist);
1800     }
1801
1802     /* ############# START The update of velocities and positions ######### */
1803     int nth = gmx_omp_nthreads_get(emntUpdate);
1804
1805 #pragma omp parallel for num_threads(nth) schedule(static)
1806     for (int th = 0; th < nth; th++)
1807     {
1808         try
1809         {
1810             int start_th, end_th;
1811             getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
1812
1813             const rvec *x_rvec  = state->x.rvec_array();
1814             rvec       *xp_rvec = upd->xp()->rvec_array();
1815             rvec       *v_rvec  = state->v.rvec_array();
1816             const rvec *f_rvec  = as_rvec_array(f.unpaddedArrayRef().data());
1817
1818             switch (inputrec->eI)
1819             {
1820                 case (eiMD):
1821                     do_update_md(start_th, end_th, step, dt,
1822                                  inputrec, md, ekind, state->box,
1823                                  x_rvec, xp_rvec, v_rvec, f_rvec,
1824                                  state->nosehoover_vxi.data(), M);
1825                     break;
1826                 case (eiSD1):
1827                     if (bDoConstr)
1828                     {
1829                         // With constraints, the SD update is done in 2 parts
1830                         doSDUpdateGeneral<SDUpdate::ForcesOnly>
1831                             (*upd->sd(),
1832                             start_th, end_th, dt,
1833                             inputrec->opts.acc, inputrec->opts.nFreeze,
1834                             md->invmass, md->ptype,
1835                             md->cFREEZE, md->cACC, nullptr,
1836                             x_rvec, xp_rvec, v_rvec, f_rvec,
1837                             step, inputrec->ld_seed, nullptr);
1838                     }
1839                     else
1840                     {
1841                         doSDUpdateGeneral<SDUpdate::Combined>
1842                             (*upd->sd(),
1843                             start_th, end_th, dt,
1844                             inputrec->opts.acc, inputrec->opts.nFreeze,
1845                             md->invmass, md->ptype,
1846                             md->cFREEZE, md->cACC, md->cTC,
1847                             x_rvec, xp_rvec, v_rvec, f_rvec,
1848                             step, inputrec->ld_seed,
1849                             DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1850                     }
1851                     break;
1852                 case (eiBD):
1853                     do_update_bd(start_th, end_th, dt,
1854                                  inputrec->opts.nFreeze, md->invmass, md->ptype,
1855                                  md->cFREEZE, md->cTC,
1856                                  x_rvec, xp_rvec, v_rvec, f_rvec,
1857                                  inputrec->bd_fric,
1858                                  upd->sd()->bd_rf.data(),
1859                                  step, inputrec->ld_seed, DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1860                     break;
1861                 case (eiVV):
1862                 case (eiVVAK):
1863                 {
1864                     gmx_bool bExtended = (inputrec->etc == etcNOSEHOOVER ||
1865                                           inputrec->epc == epcPARRINELLORAHMAN ||
1866                                           inputrec->epc == epcMTTK);
1867
1868                     /* assuming barostat coupled to group 0 */
1869                     real alpha = 1.0 + DIM/static_cast<real>(inputrec->opts.nrdf[0]);
1870                     switch (UpdatePart)
1871                     {
1872                         case etrtVELOCITY1:
1873                         case etrtVELOCITY2:
1874                             do_update_vv_vel(start_th, end_th, dt,
1875                                              inputrec->opts.acc, inputrec->opts.nFreeze,
1876                                              md->invmass, md->ptype,
1877                                              md->cFREEZE, md->cACC,
1878                                              v_rvec, f_rvec,
1879                                              bExtended, state->veta, alpha);
1880                             break;
1881                         case etrtPOSITION:
1882                             do_update_vv_pos(start_th, end_th, dt,
1883                                              inputrec->opts.nFreeze,
1884                                              md->ptype, md->cFREEZE,
1885                                              x_rvec, xp_rvec, v_rvec,
1886                                              bExtended, state->veta);
1887                             break;
1888                     }
1889                     break;
1890                 }
1891                 default:
1892                     gmx_fatal(FARGS, "Don't know how to update coordinates");
1893             }
1894         }
1895         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
1896     }
1897
1898 }
1899
1900 extern gmx_bool update_randomize_velocities(const t_inputrec *ir, int64_t step, const t_commrec *cr,
1901                                             const t_mdatoms *md,
1902                                             gmx::ArrayRef<gmx::RVec> v,
1903                                             const Update *upd,
1904                                             const gmx::Constraints *constr)
1905 {
1906
1907     real rate = (ir->delta_t)/ir->opts.tau_t[0];
1908
1909     if (ir->etc == etcANDERSEN && constr != nullptr)
1910     {
1911         /* Currently, Andersen thermostat does not support constrained
1912            systems. Functionality exists in the andersen_tcoupl
1913            function in GROMACS 4.5.7 to allow this combination. That
1914            code could be ported to the current random-number
1915            generation approach, but has not yet been done because of
1916            lack of time and resources. */
1917         gmx_fatal(FARGS, "Normal Andersen is currently not supported with constraints, use massive Andersen instead");
1918     }
1919
1920     /* proceed with andersen if 1) it's fixed probability per
1921        particle andersen or 2) it's massive andersen and it's tau_t/dt */
1922     if ((ir->etc == etcANDERSEN) || do_per_step(step, roundToInt(1.0/rate)))
1923     {
1924         andersen_tcoupl(ir, step, cr, md, v, rate,
1925                         upd->sd()->randomize_group,
1926                         upd->sd()->boltzfac);
1927         return TRUE;
1928     }
1929     return FALSE;
1930 }