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