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