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