src/gromacs/mdlib/update.cpp

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