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