src/gromacs/mdlib/update.cpp

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