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