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 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "update.h"
  41
  42 #include <cmath>
  43 #include <cstdio>
  44
  45 #include <algorithm>
  46 #include <memory>
  47
  48 #include "gromacs/domdec/domdec_struct.h"
  49 #include "gromacs/fileio/confio.h"
  50 #include "gromacs/gmxlib/network.h"
  51 #include "gromacs/gmxlib/nrnb.h"
  52 #include "gromacs/listed_forces/disre.h"
  53 #include "gromacs/listed_forces/orires.h"
  54 #include "gromacs/math/functions.h"
  55 #include "gromacs/math/invertmatrix.h"
  56 #include "gromacs/math/paddedvector.h"
  57 #include "gromacs/math/units.h"
  58 #include "gromacs/math/vec.h"
  59 #include "gromacs/math/vecdump.h"
  60 #include "gromacs/mdlib/boxdeformation.h"
  61 #include "gromacs/mdlib/constr.h"
  62 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  63 #include "gromacs/mdlib/mdatoms.h"
  64 #include "gromacs/mdlib/stat.h"
  65 #include "gromacs/mdlib/tgroup.h"
  66 #include "gromacs/mdtypes/commrec.h"
  67 #include "gromacs/mdtypes/fcdata.h"
  68 #include "gromacs/mdtypes/group.h"
  69 #include "gromacs/mdtypes/inputrec.h"
  70 #include "gromacs/mdtypes/md_enums.h"
  71 #include "gromacs/mdtypes/mdatom.h"
  72 #include "gromacs/mdtypes/state.h"
  73 #include "gromacs/pbcutil/boxutilities.h"
  74 #include "gromacs/pbcutil/pbc.h"
  75 #include "gromacs/pulling/pull.h"
  76 #include "gromacs/random/tabulatednormaldistribution.h"
  77 #include "gromacs/random/threefry.h"
  78 #include "gromacs/simd/simd.h"
  79 #include "gromacs/timing/wallcycle.h"
  80 #include "gromacs/topology/atoms.h"
  81 #include "gromacs/utility/exceptions.h"
  82 #include "gromacs/utility/fatalerror.h"
  83 #include "gromacs/utility/futil.h"
  84 #include "gromacs/utility/gmxassert.h"
  85 #include "gromacs/utility/gmxomp.h"
  86 #include "gromacs/utility/smalloc.h"
  87
  88 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
  89
  90 struct gmx_sd_const_t
  91 {
  92     double em = 0;
  93 };
  94
  95 struct gmx_sd_sigma_t
  96 {
  97     real V = 0;
  98 };
  99
 100 struct gmx_stochd_t
 101 {
 102     /* BD stuff */
 103     std::vector<real> bd_rf;
 104     /* SD stuff */
 105     std::vector<gmx_sd_const_t> sdc;
 106     std::vector<gmx_sd_sigma_t> sdsig;
 107     /* andersen temperature control stuff */
 108     std::vector<bool> randomize_group;
 109     std::vector<real> boltzfac;
 110
 111     explicit gmx_stochd_t(const t_inputrec& inputRecord);
 112 };
 113
 114 //! pImpled implementation for Update
 115 class Update::Impl
 116 {
 117 public:
 118     //! Constructor
 119     Impl(const t_inputrec& inputRecord, BoxDeformation* boxDeformation);
 120     //! Destructor
 121     ~Impl() = default;
 122
 123     void update_coords(const t_inputrec&                                inputRecord,
 124                        int64_t                                          step,
 125                        const t_mdatoms*                                 md,
 126                        t_state*                                         state,
 127                        const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
 128                        const t_fcdata&                                  fcdata,
 129                        const gmx_ekindata_t*                            ekind,
 130                        const matrix                                     M,
 131                        int                                              UpdatePart,
 132                        const t_commrec*                                 cr,
 133                        bool                                             haveConstraints);
 134
 135     void finish_update(const t_inputrec& inputRecord,
 136                        const t_mdatoms*  md,
 137                        t_state*          state,
 138                        gmx_wallcycle_t   wcycle,
 139                        bool              haveConstraints);
 140
 141     void update_sd_second_half(const t_inputrec& inputRecord,
 142                                int64_t           step,
 143                                real*             dvdlambda,
 144                                const t_mdatoms*  md,
 145                                t_state*          state,
 146                                const t_commrec*  cr,
 147                                t_nrnb*           nrnb,
 148                                gmx_wallcycle_t   wcycle,
 149                                gmx::Constraints* constr,
 150                                bool              do_log,
 151                                bool              do_ene);
 152
 153     void update_for_constraint_virial(const t_inputrec&                                inputRecord,
 154                                       const t_mdatoms&                                 md,
 155                                       const t_state&                                   state,
 156                                       const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
 157                                       const gmx_ekindata_t&                            ekind);
 158
 159     void update_temperature_constants(const t_inputrec& inputRecord);
 160
 161     const std::vector<bool>& getAndersenRandomizeGroup() const { return sd_.randomize_group; }
 162
 163     const std::vector<real>& getBoltzmanFactor() const { return sd_.boltzfac; }
 164
 165     PaddedVector<RVec>* xp() { return &xp_; }
 166
 167     BoxDeformation* deform() const { return deform_; }
 168
 169 private:
 170     //! stochastic dynamics struct
 171     gmx_stochd_t sd_;
 172     //! xprime for constraint algorithms
 173     PaddedVector<RVec> xp_;
 174     //! Box deformation handler (or nullptr if inactive).
 175     BoxDeformation* deform_ = nullptr;
 176 };
 177
 178 Update::Update(const t_inputrec& inputRecord, BoxDeformation* boxDeformation) :
 179     impl_(new Impl(inputRecord, boxDeformation)){};
 180
 181 Update::~Update(){};
 182
 183 const std::vector<bool>& Update::getAndersenRandomizeGroup() const
 184 {
 185     return impl_->getAndersenRandomizeGroup();
 186 }
 187
 188 const std::vector<real>& Update::getBoltzmanFactor() const
 189 {
 190     return impl_->getBoltzmanFactor();
 191 }
 192
 193 PaddedVector<RVec>* Update::xp()
 194 {
 195     return impl_->xp();
 196 }
 197
 198 BoxDeformation* Update::deform() const
 199 {
 200     return impl_->deform();
 201 }
 202
 203 void Update::update_coords(const t_inputrec&                                inputRecord,
 204                            int64_t                                          step,
 205                            const t_mdatoms*                                 md,
 206                            t_state*                                         state,
 207                            const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
 208                            const t_fcdata&                                  fcdata,
 209                            const gmx_ekindata_t*                            ekind,
 210                            const matrix                                     M,
 211                            int                                              updatePart,
 212                            const t_commrec*                                 cr,
 213                            const bool                                       haveConstraints)
 214 {
 215     return impl_->update_coords(inputRecord, step, md, state, f, fcdata, ekind, M, updatePart, cr,
 216                                 haveConstraints);
 217 }
 218
 219 void Update::finish_update(const t_inputrec& inputRecord,
 220                            const t_mdatoms*  md,
 221                            t_state*          state,
 222                            gmx_wallcycle_t   wcycle,
 223                            const bool        haveConstraints)
 224 {
 225     return impl_->finish_update(inputRecord, md, state, wcycle, haveConstraints);
 226 }
 227
 228 void Update::update_sd_second_half(const t_inputrec& inputRecord,
 229                                    int64_t           step,
 230                                    real*             dvdlambda,
 231                                    const t_mdatoms*  md,
 232                                    t_state*          state,
 233                                    const t_commrec*  cr,
 234                                    t_nrnb*           nrnb,
 235                                    gmx_wallcycle_t   wcycle,
 236                                    gmx::Constraints* constr,
 237                                    bool              do_log,
 238                                    bool              do_ene)
 239 {
 240     return impl_->update_sd_second_half(inputRecord, step, dvdlambda, md, state, cr, nrnb, wcycle,
 241                                         constr, do_log, do_ene);
 242 }
 243
 244 void Update::update_for_constraint_virial(const t_inputrec& inputRecord,
 245                                           const t_mdatoms&  md,
 246                                           const t_state&    state,
 247                                           const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
 248                                           const gmx_ekindata_t&                            ekind)
 249 {
 250     return impl_->update_for_constraint_virial(inputRecord, md, state, f, ekind);
 251 }
 252
 253 void Update::update_temperature_constants(const t_inputrec& inputRecord)
 254 {
 255     return impl_->update_temperature_constants(inputRecord);
 256 }
 257
 258 /*! \brief Sets the velocities of virtual sites to zero */
 259 static void clearVsiteVelocities(int start, int nrend, const unsigned short* particleType, rvec* gmx_restrict v)
 260 {
 261     for (int a = start; a < nrend; a++)
 262     {
 263         if (particleType[a] == eptVSite)
 264         {
 265             clear_rvec(v[a]);
 266         }
 267     }
 268 }
 269
 270 /*! \brief Sets whether we store the updated velocities */
 271 enum class StoreUpdatedVelocities
 272 {
 273     yes, //!< Store the updated velocities
 274     no   //!< Do not store the updated velocities
 275 };
 276
 277 /*! \brief Sets the number of different temperature coupling values */
 278 enum class NumTempScaleValues
 279 {
 280     single,  //!< Single T-scaling value (either one group or all values =1)
 281     multiple //!< Multiple T-scaling values, need to use T-group indices
 282 };
 283
 284 /*! \brief Sets if to apply no or diagonal Parrinello-Rahman pressure scaling
 285  *
 286  * Note that this enum is only used in updateMDLeapfrogSimple(), which does
 287  * not handle fully anistropic Parrinello-Rahman scaling, so we only have
 288  * options \p no and \p diagonal here and no anistropic option.
 289  */
 290 enum class ApplyParrinelloRahmanVScaling
 291 {
 292     no,      //!< Do not apply velocity scaling (not a PR-coupling run or step)
 293     diagonal //!< Apply velocity scaling using a diagonal matrix
 294 };
 295
 296 /*! \brief Integrate using leap-frog with T-scaling and optionally diagonal Parrinello-Rahman p-coupling
 297  *
 298  * \tparam       storeUpdatedVelocities Tells whether we should store the updated velocities
 299  * \tparam       numTempScaleValues     The number of different T-couple values
 300  * \tparam       applyPRVScaling        Apply Parrinello-Rahman velocity scaling
 301  * \param[in]    start                  Index of first atom to update
 302  * \param[in]    nrend                  Last atom to update: \p nrend - 1
 303  * \param[in]    dt                     The time step
 304  * \param[in]    dtPressureCouple       Time step for pressure coupling
 305  * \param[in]    invMassPerDim          1/mass per atom and dimension
 306  * \param[in]    tcstat                 Temperature coupling information
 307  * \param[in]    cTC                    T-coupling group index per atom
 308  * \param[in]    pRVScaleMatrixDiagonal Parrinello-Rahman v-scale matrix diagonal
 309  * \param[in]    x                      Input coordinates
 310  * \param[out]   xprime                 Updated coordinates
 311  * \param[inout] v                      Velocities
 312  * \param[in]    f                      Forces
 313  *
 314  * We expect this template to get good SIMD acceleration by most compilers,
 315  * unlike the more complex general template.
 316  * Note that we might get even better SIMD acceleration when we introduce
 317  * aligned (and padded) memory, possibly with some hints for the compilers.
 318  */
 319 template<StoreUpdatedVelocities storeUpdatedVelocities, NumTempScaleValues numTempScaleValues, ApplyParrinelloRahmanVScaling applyPRVScaling>
 320 static void updateMDLeapfrogSimple(int         start,
 321                                    int         nrend,
 322                                    real        dt,
 323                                    real        dtPressureCouple,
 324                                    const rvec* gmx_restrict          invMassPerDim,
 325                                    gmx::ArrayRef<const t_grp_tcstat> tcstat,
 326                                    const unsigned short*             cTC,
 327                                    const rvec                        pRVScaleMatrixDiagonal,
 328                                    const rvec* gmx_restrict x,
 329                                    rvec* gmx_restrict xprime,
 330                                    rvec* gmx_restrict v,
 331                                    const rvec* gmx_restrict f)
 332 {
 333     real lambdaGroup;
 334
 335     if (numTempScaleValues == NumTempScaleValues::single)
 336     {
 337         lambdaGroup = tcstat[0].lambda;
 338     }
 339
 340     for (int a = start; a < nrend; a++)
 341     {
 342         if (numTempScaleValues == NumTempScaleValues::multiple)
 343         {
 344             lambdaGroup = tcstat[cTC[a]].lambda;
 345         }
 346
 347         for (int d = 0; d < DIM; d++)
 348         {
 349             /* Note that using rvec invMassPerDim results in more efficient
 350              * SIMD code, but this increases the cache pressure.
 351              * For large systems with PME on the CPU this slows down the
 352              * (then already slow) update by 20%. If all data remains in cache,
 353              * using rvec is much faster.
 354              */
 355             real vNew = lambdaGroup * v[a][d] + f[a][d] * invMassPerDim[a][d] * dt;
 356
 357             if (applyPRVScaling == ApplyParrinelloRahmanVScaling::diagonal)
 358             {
 359                 vNew -= dtPressureCouple * pRVScaleMatrixDiagonal[d] * v[a][d];
 360             }
 361             if (storeUpdatedVelocities == StoreUpdatedVelocities::yes)
 362             {
 363                 v[a][d] = vNew;
 364             }
 365             xprime[a][d] = x[a][d] + vNew * dt;
 366         }
 367     }
 368 }
 369
 370 #if GMX_SIMD && GMX_SIMD_HAVE_REAL
 371 #    define GMX_HAVE_SIMD_UPDATE 1
 372 #else
 373 #    define GMX_HAVE_SIMD_UPDATE 0
 374 #endif
 375
 376 #if GMX_HAVE_SIMD_UPDATE
 377
 378 /*! \brief Load (aligned) the contents of GMX_SIMD_REAL_WIDTH rvec elements sequentially into 3 SIMD registers
 379  *
 380  * The loaded output is:
 381  * \p r0: { r[index][XX], r[index][YY], ... }
 382  * \p r1: { ... }
 383  * \p r2: { ..., r[index+GMX_SIMD_REAL_WIDTH-1][YY], r[index+GMX_SIMD_REAL_WIDTH-1][ZZ] }
 384  *
 385  * \param[in]  r      Real to an rvec array, has to be aligned to SIMD register width
 386  * \param[in]  index  Index of the first rvec triplet of reals to load
 387  * \param[out] r0     Pointer to first SIMD register
 388  * \param[out] r1     Pointer to second SIMD register
 389  * \param[out] r2     Pointer to third SIMD register
 390  */
 391 static inline void simdLoadRvecs(const rvec* r, int index, SimdReal* r0, SimdReal* r1, SimdReal* r2)
 392 {
 393     const real* realPtr = r[index];
 394
 395     GMX_ASSERT(isSimdAligned(realPtr), "Pointer should be SIMD aligned");
 396
 397     *r0 = simdLoad(realPtr + 0 * GMX_SIMD_REAL_WIDTH);
 398     *r1 = simdLoad(realPtr + 1 * GMX_SIMD_REAL_WIDTH);
 399     *r2 = simdLoad(realPtr + 2 * GMX_SIMD_REAL_WIDTH);
 400 }
 401
 402 /*! \brief Store (aligned) 3 SIMD registers sequentially to GMX_SIMD_REAL_WIDTH rvec elements
 403  *
 404  * The stored output is:
 405  * \p r[index] = { { r0[0], r0[1], ... }
 406  * ...
 407  * \p r[index+GMX_SIMD_REAL_WIDTH-1] =  { ... , r2[GMX_SIMD_REAL_WIDTH-2], r2[GMX_SIMD_REAL_WIDTH-1] }
 408  *
 409  * \param[out] r      Pointer to an rvec array, has to be aligned to SIMD register width
 410  * \param[in]  index  Index of the first rvec triplet of reals to store to
 411  * \param[in]  r0     First SIMD register
 412  * \param[in]  r1     Second SIMD register
 413  * \param[in]  r2     Third SIMD register
 414  */
 415 static inline void simdStoreRvecs(rvec* r, int index, SimdReal r0, SimdReal r1, SimdReal r2)
 416 {
 417     real* realPtr = r[index];
 418
 419     GMX_ASSERT(isSimdAligned(realPtr), "Pointer should be SIMD aligned");
 420
 421     store(realPtr + 0 * GMX_SIMD_REAL_WIDTH, r0);
 422     store(realPtr + 1 * GMX_SIMD_REAL_WIDTH, r1);
 423     store(realPtr + 2 * GMX_SIMD_REAL_WIDTH, r2);
 424 }
 425
 426 /*! \brief Integrate using leap-frog with single group T-scaling and SIMD
 427  *
 428  * \tparam       storeUpdatedVelocities Tells whether we should store the updated velocities
 429  * \param[in]    start                  Index of first atom to update
 430  * \param[in]    nrend                  Last atom to update: \p nrend - 1
 431  * \param[in]    dt                     The time step
 432  * \param[in]    invMass                1/mass per atom
 433  * \param[in]    tcstat                 Temperature coupling information
 434  * \param[in]    x                      Input coordinates
 435  * \param[out]   xprime                 Updated coordinates
 436  * \param[inout] v                      Velocities
 437  * \param[in]    f                      Forces
 438  */
 439 template<StoreUpdatedVelocities storeUpdatedVelocities>
 440 static void updateMDLeapfrogSimpleSimd(int         start,
 441                                        int         nrend,
 442                                        real        dt,
 443                                        const real* gmx_restrict          invMass,
 444                                        gmx::ArrayRef<const t_grp_tcstat> tcstat,
 445                                        const rvec* gmx_restrict x,
 446                                        rvec* gmx_restrict xprime,
 447                                        rvec* gmx_restrict v,
 448                                        const rvec* gmx_restrict f)
 449 {
 450     SimdReal timestep(dt);
 451     SimdReal lambdaSystem(tcstat[0].lambda);
 452
 453     /* We declare variables here, since code is often slower when declaring them inside the loop */
 454
 455     /* Note: We should implement a proper PaddedVector, so we don't need this check */
 456     GMX_ASSERT(isSimdAligned(invMass), "invMass should be aligned");
 457
 458     for (int a = start; a < nrend; a += GMX_SIMD_REAL_WIDTH)
 459     {
 460         SimdReal invMass0, invMass1, invMass2;
 461         expandScalarsToTriplets(simdLoad(invMass + a), &invMass0, &invMass1, &invMass2);
 462
 463         SimdReal v0, v1, v2;
 464         SimdReal f0, f1, f2;
 465         simdLoadRvecs(v, a, &v0, &v1, &v2);
 466         simdLoadRvecs(f, a, &f0, &f1, &f2);
 467
 468         v0 = fma(f0 * invMass0, timestep, lambdaSystem * v0);
 469         v1 = fma(f1 * invMass1, timestep, lambdaSystem * v1);
 470         v2 = fma(f2 * invMass2, timestep, lambdaSystem * v2);
 471
 472         if (storeUpdatedVelocities == StoreUpdatedVelocities::yes)
 473         {
 474             simdStoreRvecs(v, a, v0, v1, v2);
 475         }
 476
 477         SimdReal x0, x1, x2;
 478         simdLoadRvecs(x, a, &x0, &x1, &x2);
 479
 480         SimdReal xprime0 = fma(v0, timestep, x0);
 481         SimdReal xprime1 = fma(v1, timestep, x1);
 482         SimdReal xprime2 = fma(v2, timestep, x2);
 483
 484         simdStoreRvecs(xprime, a, xprime0, xprime1, xprime2);
 485     }
 486 }
 487
 488 #endif // GMX_HAVE_SIMD_UPDATE
 489
 490 /*! \brief Sets the NEMD acceleration type */
 491 enum class AccelerationType
 492 {
 493     none,
 494     group,
 495     cosine
 496 };
 497
 498 /*! \brief Integrate using leap-frog with support for everything.
 499  *
 500  * \tparam        accelerationType  Type of NEMD acceleration.
 501  * \param[in]     start             Index of first atom to update.
 502  * \param[in]     nrend             Last atom to update: \p nrend - 1.
 503  * \param[in]     doNoseHoover      If to apply Nose-Hoover T-coupling.
 504  * \param[in]     dt                The time step.
 505  * \param[in]     dtPressureCouple  Time step for pressure coupling, is 0 when no pressure
 506  *                                  coupling should be applied at this step.
 507  * \param[in]     accel             Acceleration per group.
 508  * \param[in]     md                Atom properties.
 509  * \param[in]     ekind             Kinetic energy data.
 510  * \param[in]     box               The box dimensions.
 511  * \param[in]     x                 Input coordinates.
 512  * \param[out]    xprime            Updated coordinates.
 513  * \param[inout]  v                 Velocities.
 514  * \param[in]     f                 Forces.
 515  * \param[in]     nh_vxi            Nose-Hoover velocity scaling factors.
 516  * \param[in]     nsttcouple        Frequency of the temperature coupling steps.
 517  * \param[in]     M                 Parrinello-Rahman scaling matrix.
 518  */
 519 template<AccelerationType accelerationType>
 520 static void updateMDLeapfrogGeneral(int                   start,
 521                                     int                   nrend,
 522                                     bool                  doNoseHoover,
 523                                     real                  dt,
 524                                     real                  dtPressureCouple,
 525                                     const rvec*           accel,
 526                                     const t_mdatoms*      md,
 527                                     const gmx_ekindata_t* ekind,
 528                                     const matrix          box,
 529                                     const rvec* gmx_restrict x,
 530                                     rvec* gmx_restrict xprime,
 531                                     rvec* gmx_restrict v,
 532                                     const rvec* gmx_restrict f,
 533                                     const double* gmx_restrict nh_vxi,
 534                                     const int                  nsttcouple,
 535                                     const matrix               M)
 536 {
 537     /* This is a version of the leap-frog integrator that supports
 538      * all combinations of T-coupling, P-coupling and NEMD.
 539      * Nose-Hoover uses the reversible leap-frog integrator from
 540      * Holian et al. Phys Rev E 52(3) : 2338, 1995
 541      */
 542
 543     gmx::ArrayRef<const t_grp_tcstat> tcstat  = ekind->tcstat;
 544     gmx::ArrayRef<const t_grp_acc>    grpstat = ekind->grpstat;
 545     const unsigned short*             cTC     = md->cTC;
 546     const unsigned short*             cACC    = md->cACC;
 547
 548     const rvec* gmx_restrict invMassPerDim = md->invMassPerDim;
 549
 550     /* Initialize group values, changed later when multiple groups are used */
 551     int ga = 0;
 552     int gt = 0;
 553
 554     real omega_Z = 2 * static_cast<real>(M_PI) / box[ZZ][ZZ];
 555
 556     for (int n = start; n < nrend; n++)
 557     {
 558         if (cTC)
 559         {
 560             gt = cTC[n];
 561         }
 562         real lg = tcstat[gt].lambda;
 563
 564         rvec vRel;
 565         real cosineZ, vCosine;
 566 #ifdef __INTEL_COMPILER
 567 #    pragma warning(disable : 280)
 568 #endif
 569         switch (accelerationType)
 570         {
 571             case AccelerationType::none: copy_rvec(v[n], vRel); break;
 572             case AccelerationType::group:
 573                 if (cACC)
 574                 {
 575                     ga = cACC[n];
 576                 }
 577                 /* Avoid scaling the group velocity */
 578                 rvec_sub(v[n], grpstat[ga].u, vRel);
 579                 break;
 580             case AccelerationType::cosine:
 581                 cosineZ = std::cos(x[n][ZZ] * omega_Z);
 582                 vCosine = cosineZ * ekind->cosacc.vcos;
 583                 /* Avoid scaling the cosine profile velocity */
 584                 copy_rvec(v[n], vRel);
 585                 vRel[XX] -= vCosine;
 586                 break;
 587         }
 588
 589         real factorNH = 0.0;
 590         if (doNoseHoover)
 591         {
 592             /* Here we account for multiple time stepping, by increasing
 593              * the Nose-Hoover correction by nsttcouple
 594              * TODO: This can be pre-computed.
 595              */
 596             factorNH = 0.5 * nsttcouple * dt * nh_vxi[gt];
 597         }
 598
 599         for (int d = 0; d < DIM; d++)
 600         {
 601             real vNew = (lg * vRel[d]
 602                          + (f[n][d] * invMassPerDim[n][d] * dt - factorNH * vRel[d]
 603                             - dtPressureCouple * iprod(M[d], vRel)))
 604                         / (1 + factorNH);
 605             switch (accelerationType)
 606             {
 607                 case AccelerationType::none: break;
 608                 case AccelerationType::group:
 609                     /* Add back the mean velocity and apply acceleration */
 610                     vNew += grpstat[ga].u[d] + accel[ga][d] * dt;
 611                     break;
 612                 case AccelerationType::cosine:
 613                     if (d == XX)
 614                     {
 615                         /* Add back the mean velocity and apply acceleration */
 616                         vNew += vCosine + cosineZ * ekind->cosacc.cos_accel * dt;
 617                     }
 618                     break;
 619             }
 620             v[n][d]      = vNew;
 621             xprime[n][d] = x[n][d] + vNew * dt;
 622         }
 623     }
 624 }
 625
 626 /*! \brief Handles the Leap-frog MD x and v integration */
 627 static void do_update_md(int         start,
 628                          int         nrend,
 629                          real        dt,
 630                          int64_t     step,
 631                          const rvec* gmx_restrict x,
 632                          rvec* gmx_restrict xprime,
 633                          rvec* gmx_restrict v,
 634                          const rvec* gmx_restrict f,
 635                          const rvec* gmx_restrict accel,
 636                          const int                etc,
 637                          const int                epc,
 638                          const int                nsttcouple,
 639                          const int                nstpcouple,
 640                          const t_mdatoms*         md,
 641                          const gmx_ekindata_t*    ekind,
 642                          const matrix             box,
 643                          const double* gmx_restrict nh_vxi,
 644                          const matrix               M)
 645 {
 646     GMX_ASSERT(nrend == start || xprime != x,
 647                "For SIMD optimization certain compilers need to have xprime != x");
 648
 649     /* Note: Berendsen pressure scaling is handled after do_update_md() */
 650     bool doTempCouple = (etc != etcNO && do_per_step(step + nsttcouple - 1, nsttcouple));
 651     bool doNoseHoover = (etc == etcNOSEHOOVER && doTempCouple);
 652     bool doParrinelloRahman =
 653             (epc == epcPARRINELLORAHMAN && do_per_step(step + nstpcouple - 1, nstpcouple));
 654     bool doPROffDiagonal = (doParrinelloRahman && (M[YY][XX] != 0 || M[ZZ][XX] != 0 || M[ZZ][YY] != 0));
 655
 656     real dtPressureCouple = (doParrinelloRahman ? nstpcouple * dt : 0);
 657
 658     /* NEMD (also cosine) acceleration is applied in updateMDLeapFrogGeneral */
 659     bool doAcceleration = (ekind->bNEMD || ekind->cosacc.cos_accel != 0);
 660
 661     if (doNoseHoover || doPROffDiagonal || doAcceleration)
 662     {
 663         matrix stepM;
 664         if (!doParrinelloRahman)
 665         {
 666             /* We should not apply PR scaling at this step */
 667             clear_mat(stepM);
 668         }
 669         else
 670         {
 671             copy_mat(M, stepM);
 672         }
 673
 674         if (!doAcceleration)
 675         {
 676             updateMDLeapfrogGeneral<AccelerationType::none>(start, nrend, doNoseHoover, dt,
 677                                                             dtPressureCouple, accel, md, ekind, box, x,
 678                                                             xprime, v, f, nh_vxi, nsttcouple, stepM);
 679         }
 680         else if (ekind->bNEMD)
 681         {
 682             updateMDLeapfrogGeneral<AccelerationType::group>(
 683                     start, nrend, doNoseHoover, dt, dtPressureCouple, accel, md, ekind, box, x,
 684                     xprime, v, f, nh_vxi, nsttcouple, stepM);
 685         }
 686         else
 687         {
 688             updateMDLeapfrogGeneral<AccelerationType::cosine>(
 689                     start, nrend, doNoseHoover, dt, dtPressureCouple, accel, md, ekind, box, x,
 690                     xprime, v, f, nh_vxi, nsttcouple, stepM);
 691         }
 692     }
 693     else
 694     {
 695         /* Use a simple and thus more efficient integration loop. */
 696         /* The simple loop does not check for particle type (so it can
 697          * be vectorized), which means we need to clear the velocities
 698          * of virtual sites in advance, when present. Note that vsite
 699          * velocities are computed after update and constraints from
 700          * their displacement.
 701          */
 702         if (md->haveVsites)
 703         {
 704             /* Note: The overhead of this loop is completely neligible */
 705             clearVsiteVelocities(start, nrend, md->ptype, v);
 706         }
 707
 708         /* We determine if we have a single T-coupling lambda value for all
 709          * atoms. That allows for better SIMD acceleration in the template.
 710          * If we do not do temperature coupling (in the run or this step),
 711          * all scaling values are 1, so we effectively have a single value.
 712          */
 713         bool haveSingleTempScaleValue = (!doTempCouple || ekind->ngtc == 1);
 714
 715         /* Extract some pointers needed by all cases */
 716         const unsigned short*             cTC           = md->cTC;
 717         gmx::ArrayRef<const t_grp_tcstat> tcstat        = ekind->tcstat;
 718         const rvec*                       invMassPerDim = md->invMassPerDim;
 719
 720         if (doParrinelloRahman)
 721         {
 722             GMX_ASSERT(!doPROffDiagonal,
 723                        "updateMDLeapfrogSimple only support diagonal Parrinello-Rahman scaling "
 724                        "matrices");
 725
 726             rvec diagM;
 727             for (int d = 0; d < DIM; d++)
 728             {
 729                 diagM[d] = M[d][d];
 730             }
 731
 732             if (haveSingleTempScaleValue)
 733             {
 734                 updateMDLeapfrogSimple<StoreUpdatedVelocities::yes, NumTempScaleValues::single,
 735                                        ApplyParrinelloRahmanVScaling::diagonal>(
 736                         start, nrend, dt, dtPressureCouple, invMassPerDim, tcstat, cTC, diagM, x,
 737                         xprime, v, f);
 738             }
 739             else
 740             {
 741                 updateMDLeapfrogSimple<StoreUpdatedVelocities::yes, NumTempScaleValues::multiple,
 742                                        ApplyParrinelloRahmanVScaling::diagonal>(
 743                         start, nrend, dt, dtPressureCouple, invMassPerDim, tcstat, cTC, diagM, x,
 744                         xprime, v, f);
 745             }
 746         }
 747         else
 748         {
 749             if (haveSingleTempScaleValue)
 750             {
 751                 /* Note that modern compilers are pretty good at vectorizing
 752                  * updateMDLeapfrogSimple(). But the SIMD version will still
 753                  * be faster because invMass lowers the cache pressure
 754                  * compared to invMassPerDim.
 755                  */
 756 #if GMX_HAVE_SIMD_UPDATE
 757                 /* Check if we can use invmass instead of invMassPerDim */
 758                 if (!md->havePartiallyFrozenAtoms)
 759                 {
 760                     updateMDLeapfrogSimpleSimd<StoreUpdatedVelocities::yes>(
 761                             start, nrend, dt, md->invmass, tcstat, x, xprime, v, f);
 762                 }
 763                 else
 764 #endif
 765                 {
 766                     updateMDLeapfrogSimple<StoreUpdatedVelocities::yes, NumTempScaleValues::single,
 767                                            ApplyParrinelloRahmanVScaling::no>(
 768                             start, nrend, dt, dtPressureCouple, invMassPerDim, tcstat, cTC, nullptr,
 769                             x, xprime, v, f);
 770                 }
 771             }
 772             else
 773             {
 774                 updateMDLeapfrogSimple<StoreUpdatedVelocities::yes, NumTempScaleValues::multiple,
 775                                        ApplyParrinelloRahmanVScaling::no>(
 776                         start, nrend, dt, dtPressureCouple, invMassPerDim, tcstat, cTC, nullptr, x,
 777                         xprime, v, f);
 778             }
 779         }
 780     }
 781 }
 782 /*! \brief Handles the Leap-frog MD x and v integration */
 783 static void doUpdateMDDoNotUpdateVelocities(int         start,
 784                                             int         nrend,
 785                                             real        dt,
 786                                             const rvec* gmx_restrict x,
 787                                             rvec* gmx_restrict xprime,
 788                                             rvec* gmx_restrict v,
 789                                             const rvec* gmx_restrict f,
 790                                             const t_mdatoms&         md,
 791                                             const gmx_ekindata_t&    ekind)
 792 {
 793     GMX_ASSERT(nrend == start || xprime != x,
 794                "For SIMD optimization certain compilers need to have xprime != x");
 795
 796     gmx::ArrayRef<const t_grp_tcstat> tcstat = ekind.tcstat;
 797
 798     /* Check if we can use invmass instead of invMassPerDim */
 799 #if GMX_HAVE_SIMD_UPDATE
 800     if (!md.havePartiallyFrozenAtoms)
 801     {
 802         updateMDLeapfrogSimpleSimd<StoreUpdatedVelocities::no>(start, nrend, dt, md.invmass, tcstat,
 803                                                                x, xprime, v, f);
 804     }
 805     else
 806 #endif
 807     {
 808         updateMDLeapfrogSimple<StoreUpdatedVelocities::no, NumTempScaleValues::single, ApplyParrinelloRahmanVScaling::no>(
 809                 start, nrend, dt, dt, md.invMassPerDim, tcstat, nullptr, nullptr, x, xprime, v, f);
 810     }
 811 }
 812
 813 static void do_update_vv_vel(int                  start,
 814                              int                  nrend,
 815                              real                 dt,
 816                              const rvec           accel[],
 817                              const ivec           nFreeze[],
 818                              const real           invmass[],
 819                              const unsigned short ptype[],
 820                              const unsigned short cFREEZE[],
 821                              const unsigned short cACC[],
 822                              rvec                 v[],
 823                              const rvec           f[],
 824                              gmx_bool             bExtended,
 825                              real                 veta,
 826                              real                 alpha)
 827 {
 828     int  gf = 0, ga = 0;
 829     int  n, d;
 830     real g, mv1, mv2;
 831
 832     if (bExtended)
 833     {
 834         g   = 0.25 * dt * veta * alpha;
 835         mv1 = std::exp(-g);
 836         mv2 = gmx::series_sinhx(g);
 837     }
 838     else
 839     {
 840         mv1 = 1.0;
 841         mv2 = 1.0;
 842     }
 843     for (n = start; n < nrend; n++)
 844     {
 845         real w_dt = invmass[n] * dt;
 846         if (cFREEZE)
 847         {
 848             gf = cFREEZE[n];
 849         }
 850         if (cACC)
 851         {
 852             ga = cACC[n];
 853         }
 854
 855         for (d = 0; d < DIM; d++)
 856         {
 857             if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
 858             {
 859                 v[n][d] = mv1 * (mv1 * v[n][d] + 0.5 * (w_dt * mv2 * f[n][d])) + 0.5 * accel[ga][d] * dt;
 860             }
 861             else
 862             {
 863                 v[n][d] = 0.0;
 864             }
 865         }
 866     }
 867 } /* do_update_vv_vel */
 868
 869 static void do_update_vv_pos(int                  start,
 870                              int                  nrend,
 871                              real                 dt,
 872                              const ivec           nFreeze[],
 873                              const unsigned short ptype[],
 874                              const unsigned short cFREEZE[],
 875                              const rvec           x[],
 876                              rvec                 xprime[],
 877                              const rvec           v[],
 878                              gmx_bool             bExtended,
 879                              real                 veta)
 880 {
 881     int  gf = 0;
 882     int  n, d;
 883     real g, mr1, mr2;
 884
 885     /* Would it make more sense if Parrinello-Rahman was put here? */
 886     if (bExtended)
 887     {
 888         g   = 0.5 * dt * veta;
 889         mr1 = std::exp(g);
 890         mr2 = gmx::series_sinhx(g);
 891     }
 892     else
 893     {
 894         mr1 = 1.0;
 895         mr2 = 1.0;
 896     }
 897
 898     for (n = start; n < nrend; n++)
 899     {
 900
 901         if (cFREEZE)
 902         {
 903             gf = cFREEZE[n];
 904         }
 905
 906         for (d = 0; d < DIM; d++)
 907         {
 908             if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
 909             {
 910                 xprime[n][d] = mr1 * (mr1 * x[n][d] + mr2 * dt * v[n][d]);
 911             }
 912             else
 913             {
 914                 xprime[n][d] = x[n][d];
 915             }
 916         }
 917     }
 918 } /* do_update_vv_pos */
 919
 920 gmx_stochd_t::gmx_stochd_t(const t_inputrec& inputRecord)
 921 {
 922     const t_grpopts* opts = &inputRecord.opts;
 923     const int        ngtc = opts->ngtc;
 924
 925     if (inputRecord.eI == eiBD)
 926     {
 927         bd_rf.resize(ngtc);
 928     }
 929     else if (EI_SD(inputRecord.eI))
 930     {
 931         sdc.resize(ngtc);
 932         sdsig.resize(ngtc);
 933
 934         for (int gt = 0; gt < ngtc; gt++)
 935         {
 936             if (opts->tau_t[gt] > 0)
 937             {
 938                 sdc[gt].em = std::exp(-inputRecord.delta_t / opts->tau_t[gt]);
 939             }
 940             else
 941             {
 942                 /* No friction and noise on this group */
 943                 sdc[gt].em = 1;
 944             }
 945         }
 946     }
 947     else if (ETC_ANDERSEN(inputRecord.etc))
 948     {
 949         randomize_group.resize(ngtc);
 950         boltzfac.resize(ngtc);
 951
 952         /* for now, assume that all groups, if randomized, are randomized at the same rate, i.e. tau_t is the same. */
 953         /* since constraint groups don't necessarily match up with temperature groups! This is checked in readir.c */
 954
 955         for (int gt = 0; gt < ngtc; gt++)
 956         {
 957             real reft = std::max<real>(0, opts->ref_t[gt]);
 958             if ((opts->tau_t[gt] > 0)
 959                 && (reft > 0)) /* tau_t or ref_t = 0 means that no randomization is done */
 960             {
 961                 randomize_group[gt] = true;
 962                 boltzfac[gt]        = BOLTZ * opts->ref_t[gt];
 963             }
 964             else
 965             {
 966                 randomize_group[gt] = false;
 967             }
 968         }
 969     }
 970 }
 971
 972 void Update::Impl::update_temperature_constants(const t_inputrec& inputRecord)
 973 {
 974     if (inputRecord.eI == eiBD)
 975     {
 976         if (inputRecord.bd_fric != 0)
 977         {
 978             for (int gt = 0; gt < inputRecord.opts.ngtc; gt++)
 979             {
 980                 sd_.bd_rf[gt] = std::sqrt(2.0 * BOLTZ * inputRecord.opts.ref_t[gt]
 981                                           / (inputRecord.bd_fric * inputRecord.delta_t));
 982             }
 983         }
 984         else
 985         {
 986             for (int gt = 0; gt < inputRecord.opts.ngtc; gt++)
 987             {
 988                 sd_.bd_rf[gt] = std::sqrt(2.0 * BOLTZ * inputRecord.opts.ref_t[gt]);
 989             }
 990         }
 991     }
 992     if (inputRecord.eI == eiSD1)
 993     {
 994         for (int gt = 0; gt < inputRecord.opts.ngtc; gt++)
 995         {
 996             real kT = BOLTZ * inputRecord.opts.ref_t[gt];
 997             /* The mass is accounted for later, since this differs per atom */
 998             sd_.sdsig[gt].V = std::sqrt(kT * (1 - sd_.sdc[gt].em * sd_.sdc[gt].em));
 999         }
1000     }
1001 }
1002
1003 Update::Impl::Impl(const t_inputrec& inputRecord, BoxDeformation* boxDeformation) :
1004     sd_(inputRecord),
1005     deform_(boxDeformation)
1006 {
1007     update_temperature_constants(inputRecord);
1008     xp_.resizeWithPadding(0);
1009 }
1010
1011 void Update::setNumAtoms(int numAtoms)
1012 {
1013
1014     impl_->xp()->resizeWithPadding(numAtoms);
1015 }
1016
1017 /*! \brief Sets the SD update type */
1018 enum class SDUpdate : int
1019 {
1020     ForcesOnly,
1021     FrictionAndNoiseOnly,
1022     Combined
1023 };
1024
1025 /*! \brief SD integrator update
1026  *
1027  * Two phases are required in the general case of a constrained
1028  * update, the first phase from the contribution of forces, before
1029  * applying constraints, and then a second phase applying the friction
1030  * and noise, and then further constraining. For details, see
1031  * Goga2012.
1032  *
1033  * Without constraints, the two phases can be combined, for
1034  * efficiency.
1035  *
1036  * Thus three instantiations of this templated function will be made,
1037  * two with only one contribution, and one with both contributions. */
1038 template<SDUpdate updateType>
1039 static void doSDUpdateGeneral(const gmx_stochd_t&  sd,
1040                               int                  start,
1041                               int                  nrend,
1042                               real                 dt,
1043                               const rvec           accel[],
1044                               const ivec           nFreeze[],
1045                               const real           invmass[],
1046                               const unsigned short ptype[],
1047                               const unsigned short cFREEZE[],
1048                               const unsigned short cACC[],
1049                               const unsigned short cTC[],
1050                               const rvec           x[],
1051                               rvec                 xprime[],
1052                               rvec                 v[],
1053                               const rvec           f[],
1054                               int64_t              step,
1055                               int                  seed,
1056                               const int*           gatindex)
1057 {
1058     // cTC, cACC and cFREEZE can be nullptr any time, but various
1059     // instantiations do not make sense with particular pointer
1060     // values.
1061     if (updateType == SDUpdate::ForcesOnly)
1062     {
1063         GMX_ASSERT(f != nullptr, "SD update with only forces requires forces");
1064         GMX_ASSERT(cTC == nullptr, "SD update with only forces cannot handle temperature groups");
1065     }
1066     if (updateType == SDUpdate::FrictionAndNoiseOnly)
1067     {
1068         GMX_ASSERT(f == nullptr, "SD update with only noise cannot handle forces");
1069         GMX_ASSERT(cACC == nullptr, "SD update with only noise cannot handle acceleration groups");
1070     }
1071     if (updateType == SDUpdate::Combined)
1072     {
1073         GMX_ASSERT(f != nullptr, "SD update with forces and noise requires forces");
1074     }
1075
1076     // Even 0 bits internal counter gives 2x64 ints (more than enough for three table lookups)
1077     gmx::ThreeFry2x64<0>                       rng(seed, gmx::RandomDomain::UpdateCoordinates);
1078     gmx::TabulatedNormalDistribution<real, 14> dist;
1079
1080     for (int n = start; n < nrend; n++)
1081     {
1082         int globalAtomIndex = gatindex ? gatindex[n] : n;
1083         rng.restart(step, globalAtomIndex);
1084         dist.reset();
1085
1086         real inverseMass = invmass[n];
1087         real invsqrtMass = std::sqrt(inverseMass);
1088
1089         int freezeGroup       = cFREEZE ? cFREEZE[n] : 0;
1090         int accelerationGroup = cACC ? cACC[n] : 0;
1091         int temperatureGroup  = cTC ? cTC[n] : 0;
1092
1093         for (int d = 0; d < DIM; d++)
1094         {
1095             if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[freezeGroup][d])
1096             {
1097                 if (updateType == SDUpdate::ForcesOnly)
1098                 {
1099                     real vn = v[n][d] + (inverseMass * f[n][d] + accel[accelerationGroup][d]) * dt;
1100                     v[n][d] = vn;
1101                     // Simple position update.
1102                     xprime[n][d] = x[n][d] + v[n][d] * dt;
1103                 }
1104                 else if (updateType == SDUpdate::FrictionAndNoiseOnly)
1105                 {
1106                     real vn = v[n][d];
1107                     v[n][d] = (vn * sd.sdc[temperatureGroup].em
1108                                + invsqrtMass * sd.sdsig[temperatureGroup].V * dist(rng));
1109                     // The previous phase already updated the
1110                     // positions with a full v*dt term that must
1111                     // now be half removed.
1112                     xprime[n][d] = xprime[n][d] + 0.5 * (v[n][d] - vn) * dt;
1113                 }
1114                 else
1115                 {
1116                     real vn = v[n][d] + (inverseMass * f[n][d] + accel[accelerationGroup][d]) * dt;
1117                     v[n][d] = (vn * sd.sdc[temperatureGroup].em
1118                                + invsqrtMass * sd.sdsig[temperatureGroup].V * dist(rng));
1119                     // Here we include half of the friction+noise
1120                     // update of v into the position update.
1121                     xprime[n][d] = x[n][d] + 0.5 * (vn + v[n][d]) * dt;
1122                 }
1123             }
1124             else
1125             {
1126                 // When using constraints, the update is split into
1127                 // two phases, but we only need to zero the update of
1128                 // virtual, shell or frozen particles in at most one
1129                 // of the phases.
1130                 if (updateType != SDUpdate::FrictionAndNoiseOnly)
1131                 {
1132                     v[n][d]      = 0.0;
1133                     xprime[n][d] = x[n][d];
1134                 }
1135             }
1136         }
1137     }
1138 }
1139
1140 static void do_update_sd(int         start,
1141                          int         nrend,
1142                          real        dt,
1143                          int64_t     step,
1144                          const rvec* gmx_restrict x,
1145                          rvec* gmx_restrict xprime,
1146                          rvec* gmx_restrict v,
1147                          const rvec* gmx_restrict f,
1148                          const rvec               accel[],
1149                          const ivec               nFreeze[],
1150                          const real               invmass[],
1151                          const unsigned short     ptype[],
1152                          const unsigned short     cFREEZE[],
1153                          const unsigned short     cACC[],
1154                          const unsigned short     cTC[],
1155                          int                      seed,
1156                          const t_commrec*         cr,
1157                          const gmx_stochd_t&      sd,
1158                          bool                     haveConstraints)
1159 {
1160     if (haveConstraints)
1161     {
1162         // With constraints, the SD update is done in 2 parts
1163         doSDUpdateGeneral<SDUpdate::ForcesOnly>(sd, start, nrend, dt, accel, nFreeze, invmass,
1164                                                 ptype, cFREEZE, cACC, nullptr, x, xprime, v, f,
1165                                                 step, seed, nullptr);
1166     }
1167     else
1168     {
1169         doSDUpdateGeneral<SDUpdate::Combined>(
1170                 sd, start, nrend, dt, accel, nFreeze, invmass, ptype, cFREEZE, cACC, cTC, x, xprime,
1171                 v, f, step, seed, DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1172     }
1173 }
1174
1175 static void do_update_bd(int         start,
1176                          int         nrend,
1177                          real        dt,
1178                          int64_t     step,
1179                          const rvec* gmx_restrict x,
1180                          rvec* gmx_restrict xprime,
1181                          rvec* gmx_restrict v,
1182                          const rvec* gmx_restrict f,
1183                          const ivec               nFreeze[],
1184                          const real               invmass[],
1185                          const unsigned short     ptype[],
1186                          const unsigned short     cFREEZE[],
1187                          const unsigned short     cTC[],
1188                          real                     friction_coefficient,
1189                          const real*              rf,
1190                          int                      seed,
1191                          const int*               gatindex)
1192 {
1193     /* note -- these appear to be full step velocities . . .  */
1194     int  gf = 0, gt = 0;
1195     real vn;
1196     real invfr = 0;
1197     int  n, d;
1198     // Use 1 bit of internal counters to give us 2*2 64-bits values per stream
1199     // Each 64-bit value is enough for 4 normal distribution table numbers.
1200     gmx::ThreeFry2x64<0>                       rng(seed, gmx::RandomDomain::UpdateCoordinates);
1201     gmx::TabulatedNormalDistribution<real, 14> dist;
1202
1203     if (friction_coefficient != 0)
1204     {
1205         invfr = 1.0 / friction_coefficient;
1206     }
1207
1208     for (n = start; (n < nrend); n++)
1209     {
1210         int ng = gatindex ? gatindex[n] : n;
1211
1212         rng.restart(step, ng);
1213         dist.reset();
1214
1215         if (cFREEZE)
1216         {
1217             gf = cFREEZE[n];
1218         }
1219         if (cTC)
1220         {
1221             gt = cTC[n];
1222         }
1223         for (d = 0; (d < DIM); d++)
1224         {
1225             if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
1226             {
1227                 if (friction_coefficient != 0)
1228                 {
1229                     vn = invfr * f[n][d] + rf[gt] * dist(rng);
1230                 }
1231                 else
1232                 {
1233                     /* NOTE: invmass = 2/(mass*friction_constant*dt) */
1234                     vn = 0.5 * invmass[n] * f[n][d] * dt
1235                          + std::sqrt(0.5 * invmass[n]) * rf[gt] * dist(rng);
1236                 }
1237
1238                 v[n][d]      = vn;
1239                 xprime[n][d] = x[n][d] + vn * dt;
1240             }
1241             else
1242             {
1243                 v[n][d]      = 0.0;
1244                 xprime[n][d] = x[n][d];
1245             }
1246         }
1247     }
1248 }
1249
1250 extern void init_ekinstate(ekinstate_t* ekinstate, const t_inputrec* ir)
1251 {
1252     ekinstate->ekin_n = ir->opts.ngtc;
1253     snew(ekinstate->ekinh, ekinstate->ekin_n);
1254     snew(ekinstate->ekinf, ekinstate->ekin_n);
1255     snew(ekinstate->ekinh_old, ekinstate->ekin_n);
1256     ekinstate->ekinscalef_nhc.resize(ekinstate->ekin_n);
1257     ekinstate->ekinscaleh_nhc.resize(ekinstate->ekin_n);
1258     ekinstate->vscale_nhc.resize(ekinstate->ekin_n);
1259     ekinstate->dekindl          = 0;
1260     ekinstate->mvcos            = 0;
1261     ekinstate->hasReadEkinState = false;
1262 }
1263
1264 void update_ekinstate(ekinstate_t* ekinstate, const gmx_ekindata_t* ekind)
1265 {
1266     int i;
1267
1268     for (i = 0; i < ekinstate->ekin_n; i++)
1269     {
1270         copy_mat(ekind->tcstat[i].ekinh, ekinstate->ekinh[i]);
1271         copy_mat(ekind->tcstat[i].ekinf, ekinstate->ekinf[i]);
1272         copy_mat(ekind->tcstat[i].ekinh_old, ekinstate->ekinh_old[i]);
1273         ekinstate->ekinscalef_nhc[i] = ekind->tcstat[i].ekinscalef_nhc;
1274         ekinstate->ekinscaleh_nhc[i] = ekind->tcstat[i].ekinscaleh_nhc;
1275         ekinstate->vscale_nhc[i]     = ekind->tcstat[i].vscale_nhc;
1276     }
1277
1278     copy_mat(ekind->ekin, ekinstate->ekin_total);
1279     ekinstate->dekindl = ekind->dekindl;
1280     ekinstate->mvcos   = ekind->cosacc.mvcos;
1281 }
1282
1283 void restore_ekinstate_from_state(const t_commrec* cr, gmx_ekindata_t* ekind, const ekinstate_t* ekinstate)
1284 {
1285     int i, n;
1286
1287     if (MASTER(cr))
1288     {
1289         for (i = 0; i < ekinstate->ekin_n; i++)
1290         {
1291             copy_mat(ekinstate->ekinh[i], ekind->tcstat[i].ekinh);
1292             copy_mat(ekinstate->ekinf[i], ekind->tcstat[i].ekinf);
1293             copy_mat(ekinstate->ekinh_old[i], ekind->tcstat[i].ekinh_old);
1294             ekind->tcstat[i].ekinscalef_nhc = ekinstate->ekinscalef_nhc[i];
1295             ekind->tcstat[i].ekinscaleh_nhc = ekinstate->ekinscaleh_nhc[i];
1296             ekind->tcstat[i].vscale_nhc     = ekinstate->vscale_nhc[i];
1297         }
1298
1299         copy_mat(ekinstate->ekin_total, ekind->ekin);
1300
1301         ekind->dekindl      = ekinstate->dekindl;
1302         ekind->cosacc.mvcos = ekinstate->mvcos;
1303         n                   = ekinstate->ekin_n;
1304     }
1305
1306     if (PAR(cr))
1307     {
1308         gmx_bcast(sizeof(n), &n, cr->mpi_comm_mygroup);
1309         for (i = 0; i < n; i++)
1310         {
1311             gmx_bcast(DIM * DIM * sizeof(ekind->tcstat[i].ekinh[0][0]), ekind->tcstat[i].ekinh[0],
1312                       cr->mpi_comm_mygroup);
1313             gmx_bcast(DIM * DIM * sizeof(ekind->tcstat[i].ekinf[0][0]), ekind->tcstat[i].ekinf[0],
1314                       cr->mpi_comm_mygroup);
1315             gmx_bcast(DIM * DIM * sizeof(ekind->tcstat[i].ekinh_old[0][0]),
1316                       ekind->tcstat[i].ekinh_old[0], cr->mpi_comm_mygroup);
1317
1318             gmx_bcast(sizeof(ekind->tcstat[i].ekinscalef_nhc), &(ekind->tcstat[i].ekinscalef_nhc),
1319                       cr->mpi_comm_mygroup);
1320             gmx_bcast(sizeof(ekind->tcstat[i].ekinscaleh_nhc), &(ekind->tcstat[i].ekinscaleh_nhc),
1321                       cr->mpi_comm_mygroup);
1322             gmx_bcast(sizeof(ekind->tcstat[i].vscale_nhc), &(ekind->tcstat[i].vscale_nhc),
1323                       cr->mpi_comm_mygroup);
1324         }
1325         gmx_bcast(DIM * DIM * sizeof(ekind->ekin[0][0]), ekind->ekin[0], cr->mpi_comm_mygroup);
1326
1327         gmx_bcast(sizeof(ekind->dekindl), &ekind->dekindl, cr->mpi_comm_mygroup);
1328         gmx_bcast(sizeof(ekind->cosacc.mvcos), &ekind->cosacc.mvcos, cr->mpi_comm_mygroup);
1329     }
1330 }
1331
1332 void getThreadAtomRange(int numThreads, int threadIndex, int numAtoms, int* startAtom, int* endAtom)
1333 {
1334 #if GMX_HAVE_SIMD_UPDATE
1335     constexpr int blockSize = GMX_SIMD_REAL_WIDTH;
1336 #else
1337     constexpr int blockSize = 1;
1338 #endif
1339     int numBlocks = (numAtoms + blockSize - 1) / blockSize;
1340
1341     *startAtom = ((numBlocks * threadIndex) / numThreads) * blockSize;
1342     *endAtom   = ((numBlocks * (threadIndex + 1)) / numThreads) * blockSize;
1343     if (threadIndex == numThreads - 1)
1344     {
1345         *endAtom = numAtoms;
1346     }
1347 }
1348
1349 void Update::Impl::update_sd_second_half(const t_inputrec& inputRecord,
1350                                          int64_t           step,
1351                                          real*             dvdlambda,
1352                                          const t_mdatoms*  md,
1353                                          t_state*          state,
1354                                          const t_commrec*  cr,
1355                                          t_nrnb*           nrnb,
1356                                          gmx_wallcycle_t   wcycle,
1357                                          gmx::Constraints* constr,
1358                                          bool              do_log,
1359                                          bool              do_ene)
1360 {
1361     if (!constr)
1362     {
1363         return;
1364     }
1365     if (inputRecord.eI == eiSD1)
1366     {
1367         int homenr = md->homenr;
1368
1369         /* Cast delta_t from double to real to make the integrators faster.
1370          * The only reason for having delta_t double is to get accurate values
1371          * for t=delta_t*step when step is larger than float precision.
1372          * For integration dt the accuracy of real suffices, since with
1373          * integral += dt*integrand the increment is nearly always (much) smaller
1374          * than the integral (and the integrand has real precision).
1375          */
1376         real dt = inputRecord.delta_t;
1377
1378         wallcycle_start(wcycle, ewcUPDATE);
1379
1380         int nth = gmx_omp_nthreads_get(emntUpdate);
1381
1382 #pragma omp parallel for num_threads(nth) schedule(static)
1383         for (int th = 0; th < nth; th++)
1384         {
1385             try
1386             {
1387                 int start_th, end_th;
1388                 getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
1389
1390                 doSDUpdateGeneral<SDUpdate::FrictionAndNoiseOnly>(
1391                         sd_, start_th, end_th, dt, inputRecord.opts.acc, inputRecord.opts.nFreeze,
1392                         md->invmass, md->ptype, md->cFREEZE, nullptr, md->cTC, state->x.rvec_array(),
1393                         xp_.rvec_array(), state->v.rvec_array(), nullptr, step, inputRecord.ld_seed,
1394                         DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1395             }
1396             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1397         }
1398         inc_nrnb(nrnb, eNR_UPDATE, homenr);
1399         wallcycle_stop(wcycle, ewcUPDATE);
1400
1401         /* Constrain the coordinates upd->xp for half a time step */
1402         bool computeVirial = false;
1403         constr->apply(do_log, do_ene, step, 1, 0.5, state->x.arrayRefWithPadding(),
1404                       xp_.arrayRefWithPadding(), ArrayRef<RVec>(), state->box,
1405                       state->lambda[efptBONDED], dvdlambda, state->v.arrayRefWithPadding(),
1406                       computeVirial, nullptr, ConstraintVariable::Positions);
1407     }
1408 }
1409
1410 void Update::Impl::finish_update(const t_inputrec& inputRecord,
1411                                  const t_mdatoms*  md,
1412                                  t_state*          state,
1413                                  gmx_wallcycle_t   wcycle,
1414                                  const bool        haveConstraints)
1415 {
1416     /* NOTE: Currently we always integrate to a temporary buffer and
1417      * then copy the results back here.
1418      */
1419
1420     wallcycle_start_nocount(wcycle, ewcUPDATE);
1421
1422     const int homenr = md->homenr;
1423     auto      xp     = makeConstArrayRef(xp_).subArray(0, homenr);
1424     auto      x      = makeArrayRef(state->x).subArray(0, homenr);
1425
1426     if (md->havePartiallyFrozenAtoms && haveConstraints)
1427     {
1428         /* We have atoms that are frozen along some, but not all dimensions,
1429          * then constraints will have moved them also along the frozen dimensions.
1430          * To freeze such degrees of freedom we do not copy them back here.
1431          */
1432         const ivec* nFreeze = inputRecord.opts.nFreeze;
1433
1434         for (int i = 0; i < homenr; i++)
1435         {
1436             const int g = md->cFREEZE[i];
1437
1438             for (int d = 0; d < DIM; d++)
1439             {
1440                 if (nFreeze[g][d] == 0)
1441                 {
1442                     x[i][d] = xp[i][d];
1443                 }
1444             }
1445         }
1446     }
1447     else
1448     {
1449         /* We have no frozen atoms or fully frozen atoms which have not
1450          * been moved by the update, so we can simply copy all coordinates.
1451          */
1452         int gmx_unused nth = gmx_omp_nthreads_get(emntUpdate);
1453 #pragma omp parallel for num_threads(nth) schedule(static)
1454         for (int i = 0; i < homenr; i++)
1455         {
1456             // Trivial statement, does not throw
1457             x[i] = xp[i];
1458         }
1459     }
1460
1461     wallcycle_stop(wcycle, ewcUPDATE);
1462 }
1463
1464 void Update::Impl::update_coords(const t_inputrec&                                inputRecord,
1465                                  int64_t                                          step,
1466                                  const t_mdatoms*                                 md,
1467                                  t_state*                                         state,
1468                                  const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
1469                                  const t_fcdata&                                  fcdata,
1470                                  const gmx_ekindata_t*                            ekind,
1471                                  const matrix                                     M,
1472                                  int                                              updatePart,
1473                                  const t_commrec*                                 cr,
1474                                  const bool                                       haveConstraints)
1475 {
1476     /* Running the velocity half does nothing except for velocity verlet */
1477     if ((updatePart == etrtVELOCITY1 || updatePart == etrtVELOCITY2) && !EI_VV(inputRecord.eI))
1478     {
1479         gmx_incons("update_coords called for velocity without VV integrator");
1480     }
1481
1482     int homenr = md->homenr;
1483
1484     /* Cast to real for faster code, no loss in precision (see comment above) */
1485     real dt = inputRecord.delta_t;
1486
1487     /* We need to update the NMR restraint history when time averaging is used */
1488     if (state->flags & (1 << estDISRE_RM3TAV))
1489     {
1490         update_disres_history(*fcdata.disres, &state->hist);
1491     }
1492     if (state->flags & (1 << estORIRE_DTAV))
1493     {
1494         update_orires_history(*fcdata.orires, &state->hist);
1495     }
1496
1497     /* ############# START The update of velocities and positions ######### */
1498     int nth = gmx_omp_nthreads_get(emntUpdate);
1499
1500 #pragma omp parallel for num_threads(nth) schedule(static)
1501     for (int th = 0; th < nth; th++)
1502     {
1503         try
1504         {
1505             int start_th, end_th;
1506             getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
1507
1508             const rvec* x_rvec  = state->x.rvec_array();
1509             rvec*       xp_rvec = xp_.rvec_array();
1510             rvec*       v_rvec  = state->v.rvec_array();
1511             const rvec* f_rvec  = as_rvec_array(f.unpaddedConstArrayRef().data());
1512
1513             switch (inputRecord.eI)
1514             {
1515                 case (eiMD):
1516                     do_update_md(start_th, end_th, dt, step, x_rvec, xp_rvec, v_rvec, f_rvec,
1517                                  inputRecord.opts.acc, inputRecord.etc, inputRecord.epc,
1518                                  inputRecord.nsttcouple, inputRecord.nstpcouple, md, ekind,
1519                                  state->box, state->nosehoover_vxi.data(), M);
1520                     break;
1521                 case (eiSD1):
1522                     do_update_sd(start_th, end_th, dt, step, x_rvec, xp_rvec, v_rvec, f_rvec,
1523                                  inputRecord.opts.acc, inputRecord.opts.nFreeze, md->invmass,
1524                                  md->ptype, md->cFREEZE, md->cACC, md->cTC, inputRecord.ld_seed, cr,
1525                                  sd_, haveConstraints);
1526                     break;
1527                 case (eiBD):
1528                     do_update_bd(start_th, end_th, dt, step, x_rvec, xp_rvec, v_rvec, f_rvec,
1529                                  inputRecord.opts.nFreeze, md->invmass, md->ptype, md->cFREEZE,
1530                                  md->cTC, inputRecord.bd_fric, sd_.bd_rf.data(), inputRecord.ld_seed,
1531                                  DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1532                     break;
1533                 case (eiVV):
1534                 case (eiVVAK):
1535                 {
1536                     gmx_bool bExtended = (inputRecord.etc == etcNOSEHOOVER || inputRecord.epc == epcPARRINELLORAHMAN
1537                                           || inputRecord.epc == epcMTTK);
1538
1539                     /* assuming barostat coupled to group 0 */
1540                     real alpha = 1.0 + DIM / static_cast<real>(inputRecord.opts.nrdf[0]);
1541                     switch (updatePart)
1542                     {
1543                         case etrtVELOCITY1:
1544                         case etrtVELOCITY2:
1545                             do_update_vv_vel(start_th, end_th, dt, inputRecord.opts.acc,
1546                                              inputRecord.opts.nFreeze, md->invmass, md->ptype, md->cFREEZE,
1547                                              md->cACC, v_rvec, f_rvec, bExtended, state->veta, alpha);
1548                             break;
1549                         case etrtPOSITION:
1550                             do_update_vv_pos(start_th, end_th, dt, inputRecord.opts.nFreeze,
1551                                              md->ptype, md->cFREEZE, x_rvec, xp_rvec, v_rvec,
1552                                              bExtended, state->veta);
1553                             break;
1554                     }
1555                     break;
1556                 }
1557                 default: gmx_fatal(FARGS, "Don't know how to update coordinates");
1558             }
1559         }
1560         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1561     }
1562 }
1563
1564 void Update::Impl::update_for_constraint_virial(const t_inputrec& inputRecord,
1565                                                 const t_mdatoms&  md,
1566                                                 const t_state&    state,
1567                                                 const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
1568                                                 const gmx_ekindata_t& ekind)
1569 {
1570     GMX_ASSERT(inputRecord.eI == eiMD || inputRecord.eI == eiSD1,
1571                "Only leap-frog is supported here");
1572
1573     // Cast to real for faster code, no loss in precision
1574     const real dt = inputRecord.delta_t;
1575
1576     const int nth = gmx_omp_nthreads_get(emntUpdate);
1577
1578 #pragma omp parallel for num_threads(nth) schedule(static)
1579     for (int th = 0; th < nth; th++)
1580     {
1581         try
1582         {
1583             int start_th, end_th;
1584             getThreadAtomRange(nth, th, md.homenr, &start_th, &end_th);
1585
1586             const rvec* x_rvec  = state.x.rvec_array();
1587             rvec*       xp_rvec = xp_.rvec_array();
1588             rvec*       v_rvec  = const_cast<rvec*>(state.v.rvec_array());
1589             const rvec* f_rvec  = as_rvec_array(f.unpaddedConstArrayRef().data());
1590
1591             doUpdateMDDoNotUpdateVelocities(start_th, end_th, dt, x_rvec, xp_rvec, v_rvec, f_rvec,
1592                                             md, ekind);
1593         }
1594         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1595     }
1596 }