src/gromacs/mdlib/update.cpp

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