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,2021, 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/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/paddedvector.h"
  55 #include "gromacs/math/units.h"
  56 #include "gromacs/math/vec.h"
  57 #include "gromacs/mdlib/boxdeformation.h"
  58 #include "gromacs/mdlib/constr.h"
  59 #include "gromacs/mdlib/gmx_omp_nthreads.h"
  60 #include "gromacs/mdlib/stat.h"
  61 #include "gromacs/mdlib/tgroup.h"
  62 #include "gromacs/mdtypes/commrec.h"
  63 #include "gromacs/mdtypes/fcdata.h"
  64 #include "gromacs/mdtypes/group.h"
  65 #include "gromacs/mdtypes/inputrec.h"
  66 #include "gromacs/mdtypes/md_enums.h"
  67 #include "gromacs/mdtypes/state.h"
  68 #include "gromacs/pbcutil/pbc.h"
  69 #include "gromacs/pulling/pull.h"
  70 #include "gromacs/random/tabulatednormaldistribution.h"
  71 #include "gromacs/random/threefry.h"
  72 #include "gromacs/simd/simd.h"
  73 #include "gromacs/timing/wallcycle.h"
  74 #include "gromacs/topology/atoms.h"
  75 #include "gromacs/utility/exceptions.h"
  76 #include "gromacs/utility/fatalerror.h"
  77 #include "gromacs/utility/futil.h"
  78 #include "gromacs/utility/gmxassert.h"
  79 #include "gromacs/utility/smalloc.h"
  80
  81 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
  82
  83 struct gmx_sd_const_t
  84 {
  85     double em = 0;
  86 };
  87
  88 struct gmx_sd_sigma_t
  89 {
  90     real V = 0;
  91 };
  92
  93 struct gmx_stochd_t
  94 {
  95     /* BD stuff */
  96     std::vector<real> bd_rf;
  97     /* SD stuff */
  98     std::vector<gmx_sd_const_t> sdc;
  99     std::vector<gmx_sd_sigma_t> sdsig;
 100     /* andersen temperature control stuff */
 101     std::vector<bool> randomize_group;
 102     std::vector<real> boltzfac;
 103
 104     explicit gmx_stochd_t(const t_inputrec& inputRecord);
 105 };
 106
 107 //! pImpled implementation for Update
 108 class Update::Impl
 109 {
 110 public:
 111     //! Constructor
 112     Impl(const t_inputrec& inputRecord, BoxDeformation* boxDeformation);
 113     //! Destructor
 114     ~Impl() = default;
 115
 116     void update_coords(const t_inputrec&                                inputRecord,
 117                        int64_t                                          step,
 118                        int                                              homenr,
 119                        bool                                             havePartiallyFrozenAtoms,
 120                        gmx::ArrayRef<const ParticleType>                ptype,
 121                        gmx::ArrayRef<const real>                        invMass,
 122                        gmx::ArrayRef<const rvec>                        invMassPerDim,
 123                        t_state*                                         state,
 124                        const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
 125                        t_fcdata*                                        fcdata,
 126                        const gmx_ekindata_t*                            ekind,
 127                        const matrix                                     M,
 128                        int                                              UpdatePart,
 129                        const t_commrec*                                 cr,
 130                        bool                                             haveConstraints);
 131
 132     void finish_update(const t_inputrec&                   inputRecord,
 133                        bool                                havePartiallyFrozenAtoms,
 134                        int                                 homenr,
 135                        gmx::ArrayRef<const unsigned short> cFREEZE,
 136                        t_state*                            state,
 137                        gmx_wallcycle*                      wcycle,
 138                        bool                                haveConstraints);
 139
 140     void update_sd_second_half(const t_inputrec&                 inputRecord,
 141                                int64_t                           step,
 142                                real*                             dvdlambda,
 143                                int                               homenr,
 144                                gmx::ArrayRef<const ParticleType> ptype,
 145                                gmx::ArrayRef<const real>         invMass,
 146                                t_state*                          state,
 147                                const t_commrec*                  cr,
 148                                t_nrnb*                           nrnb,
 149                                gmx_wallcycle*                    wcycle,
 150                                gmx::Constraints*                 constr,
 151                                bool                              do_log,
 152                                bool                              do_ene);
 153
 154     void update_for_constraint_virial(const t_inputrec&         inputRecord,
 155                                       int                       homenr,
 156                                       bool                      havePartiallyFrozenAtoms,
 157                                       gmx::ArrayRef<const real> invmass,
 158                                       gmx::ArrayRef<const rvec> invMassPerDim,
 159                                       const t_state&            state,
 160                                       const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
 161                                       const gmx_ekindata_t&                            ekind);
 162
 163     void update_temperature_constants(const t_inputrec& inputRecord);
 164
 165     const std::vector<bool>& getAndersenRandomizeGroup() const { return sd_.randomize_group; }
 166
 167     const std::vector<real>& getBoltzmanFactor() const { return sd_.boltzfac; }
 168
 169     PaddedVector<RVec>* xp() { return &xp_; }
 170
 171     BoxDeformation* deform() const { return deform_; }
 172
 173     //! Group index for freezing
 174     gmx::ArrayRef<const unsigned short> cFREEZE_;
 175     //! Group index for temperature coupling
 176     gmx::ArrayRef<const unsigned short> cTC_;
 177
 178 private:
 179     //! stochastic dynamics struct
 180     gmx_stochd_t sd_;
 181     //! xprime for constraint algorithms
 182     PaddedVector<RVec> xp_;
 183     //! Box deformation handler (or nullptr if inactive).
 184     BoxDeformation* deform_ = nullptr;
 185 };
 186
 187 Update::Update(const t_inputrec& inputRecord, BoxDeformation* boxDeformation) :
 188     impl_(new Impl(inputRecord, boxDeformation)){};
 189
 190 Update::~Update(){};
 191
 192 const std::vector<bool>& Update::getAndersenRandomizeGroup() const
 193 {
 194     return impl_->getAndersenRandomizeGroup();
 195 }
 196
 197 const std::vector<real>& Update::getBoltzmanFactor() const
 198 {
 199     return impl_->getBoltzmanFactor();
 200 }
 201
 202 PaddedVector<RVec>* Update::xp()
 203 {
 204     return impl_->xp();
 205 }
 206
 207 BoxDeformation* Update::deform() const
 208 {
 209     return impl_->deform();
 210 }
 211
 212 void Update::update_coords(const t_inputrec&                 inputRecord,
 213                            int64_t                           step,
 214                            const int                         homenr,
 215                            const bool                        havePartiallyFrozenAtoms,
 216                            gmx::ArrayRef<const ParticleType> ptype,
 217                            gmx::ArrayRef<const real>         invMass,
 218                            gmx::ArrayRef<const rvec>         invMassPerDim,
 219                            t_state*                          state,
 220                            const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
 221                            t_fcdata*                                        fcdata,
 222                            const gmx_ekindata_t*                            ekind,
 223                            const matrix                                     M,
 224                            int                                              updatePart,
 225                            const t_commrec*                                 cr,
 226                            const bool                                       haveConstraints)
 227 {
 228     return impl_->update_coords(inputRecord,
 229                                 step,
 230                                 homenr,
 231                                 havePartiallyFrozenAtoms,
 232                                 ptype,
 233                                 invMass,
 234                                 invMassPerDim,
 235                                 state,
 236                                 f,
 237                                 fcdata,
 238                                 ekind,
 239                                 M,
 240                                 updatePart,
 241                                 cr,
 242                                 haveConstraints);
 243 }
 244
 245 void Update::finish_update(const t_inputrec& inputRecord,
 246                            const bool        havePartiallyFrozenAtoms,
 247                            const int         homenr,
 248                            t_state*          state,
 249                            gmx_wallcycle*    wcycle,
 250                            const bool        haveConstraints)
 251 {
 252     return impl_->finish_update(
 253             inputRecord, havePartiallyFrozenAtoms, homenr, impl_->cFREEZE_, state, wcycle, haveConstraints);
 254 }
 255
 256 void Update::update_sd_second_half(const t_inputrec&                 inputRecord,
 257                                    int64_t                           step,
 258                                    real*                             dvdlambda,
 259                                    const int                         homenr,
 260                                    gmx::ArrayRef<const ParticleType> ptype,
 261                                    gmx::ArrayRef<const real>         invMass,
 262                                    t_state*                          state,
 263                                    const t_commrec*                  cr,
 264                                    t_nrnb*                           nrnb,
 265                                    gmx_wallcycle*                    wcycle,
 266                                    gmx::Constraints*                 constr,
 267                                    bool                              do_log,
 268                                    bool                              do_ene)
 269 {
 270     return impl_->update_sd_second_half(
 271             inputRecord, step, dvdlambda, homenr, ptype, invMass, state, cr, nrnb, wcycle, constr, do_log, do_ene);
 272 }
 273
 274 void Update::update_for_constraint_virial(const t_inputrec&         inputRecord,
 275                                           const int                 homenr,
 276                                           const bool                havePartiallyFrozenAtoms,
 277                                           gmx::ArrayRef<const real> invmass,
 278                                           gmx::ArrayRef<const rvec> invMassPerDim,
 279                                           const t_state&            state,
 280                                           const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
 281                                           const gmx_ekindata_t&                            ekind)
 282 {
 283     return impl_->update_for_constraint_virial(
 284             inputRecord, homenr, havePartiallyFrozenAtoms, invmass, invMassPerDim, state, f, ekind);
 285 }
 286
 287 void Update::update_temperature_constants(const t_inputrec& inputRecord)
 288 {
 289     return impl_->update_temperature_constants(inputRecord);
 290 }
 291
 292 /*! \brief Sets whether we store the updated velocities */
 293 enum class StoreUpdatedVelocities
 294 {
 295     yes, //!< Store the updated velocities
 296     no   //!< Do not store the updated velocities
 297 };
 298
 299 /*! \brief Sets the number of different temperature coupling values */
 300 enum class NumTempScaleValues
 301 {
 302     single,  //!< Single T-scaling value (either one group or all values =1)
 303     multiple //!< Multiple T-scaling values, need to use T-group indices
 304 };
 305
 306 /*! \brief Sets if to apply no or diagonal Parrinello-Rahman pressure scaling
 307  *
 308  * Note that this enum is only used in updateMDLeapfrogSimple(), which does
 309  * not handle fully anistropic Parrinello-Rahman scaling, so we only have
 310  * options \p no and \p diagonal here and no anistropic option.
 311  */
 312 enum class ApplyParrinelloRahmanVScaling
 313 {
 314     no,      //!< Do not apply velocity scaling (not a PR-coupling run or step)
 315     diagonal //!< Apply velocity scaling using a diagonal matrix
 316 };
 317
 318 /*! \brief Integrate using leap-frog with T-scaling and optionally diagonal Parrinello-Rahman p-coupling
 319  *
 320  * \tparam       storeUpdatedVelocities Tells whether we should store the updated velocities
 321  * \tparam       numTempScaleValues     The number of different T-couple values
 322  * \tparam       applyPRVScaling        Apply Parrinello-Rahman velocity scaling
 323  * \param[in]    start                  Index of first atom to update
 324  * \param[in]    nrend                  Last atom to update: \p nrend - 1
 325  * \param[in]    dt                     The time step
 326  * \param[in]    dtPressureCouple       Time step for pressure coupling
 327  * \param[in]    invMassPerDim          1/mass per atom and dimension
 328  * \param[in]    tcstat                 Temperature coupling information
 329  * \param[in]    cTC                    T-coupling group index per atom
 330  * \param[in]    pRVScaleMatrixDiagonal Parrinello-Rahman v-scale matrix diagonal
 331  * \param[in]    x                      Input coordinates
 332  * \param[out]   xprime                 Updated coordinates
 333  * \param[inout] v                      Velocities, type either rvec* or const rvec*
 334  * \param[in]    f                      Forces
 335  *
 336  * We expect this template to get good SIMD acceleration by most compilers,
 337  * unlike the more complex general template.
 338  * Note that we might get even better SIMD acceleration when we introduce
 339  * aligned (and padded) memory, possibly with some hints for the compilers.
 340  */
 341 template<StoreUpdatedVelocities storeUpdatedVelocities, NumTempScaleValues numTempScaleValues, ApplyParrinelloRahmanVScaling applyPRVScaling, typename VelocityType>
 342 static std::enable_if_t<std::is_same<VelocityType, rvec>::value || std::is_same<VelocityType, const rvec>::value, void>
 343 updateMDLeapfrogSimple(int                                 start,
 344                        int                                 nrend,
 345                        real                                dt,
 346                        real                                dtPressureCouple,
 347                        gmx::ArrayRef<const rvec>           invMassPerDim,
 348                        gmx::ArrayRef<const t_grp_tcstat>   tcstat,
 349                        gmx::ArrayRef<const unsigned short> cTC,
 350                        const rvec                          pRVScaleMatrixDiagonal,
 351                        const rvec* gmx_restrict            x,
 352                        rvec* gmx_restrict                  xprime,
 353                        VelocityType* gmx_restrict          v,
 354                        const rvec* gmx_restrict            f)
 355 {
 356     real lambdaGroup;
 357
 358     if (numTempScaleValues == NumTempScaleValues::single)
 359     {
 360         lambdaGroup = tcstat[0].lambda;
 361     }
 362
 363     for (int a = start; a < nrend; a++)
 364     {
 365         if (numTempScaleValues == NumTempScaleValues::multiple)
 366         {
 367             lambdaGroup = tcstat[cTC[a]].lambda;
 368         }
 369
 370         for (int d = 0; d < DIM; d++)
 371         {
 372             /* Note that using rvec invMassPerDim results in more efficient
 373              * SIMD code, but this increases the cache pressure.
 374              * For large systems with PME on the CPU this slows down the
 375              * (then already slow) update by 20%. If all data remains in cache,
 376              * using rvec is much faster.
 377              */
 378             real vNew = lambdaGroup * v[a][d] + f[a][d] * invMassPerDim[a][d] * dt;
 379
 380             if (applyPRVScaling == ApplyParrinelloRahmanVScaling::diagonal)
 381             {
 382                 vNew -= dtPressureCouple * pRVScaleMatrixDiagonal[d] * v[a][d];
 383             }
 384             if constexpr (storeUpdatedVelocities == StoreUpdatedVelocities::yes)
 385             {
 386                 v[a][d] = vNew;
 387             }
 388             // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
 389             xprime[a][d] = x[a][d] + vNew * dt;
 390         }
 391     }
 392 }
 393
 394 #if GMX_SIMD && GMX_SIMD_HAVE_REAL
 395 #    define GMX_HAVE_SIMD_UPDATE 1
 396 #else
 397 #    define GMX_HAVE_SIMD_UPDATE 0
 398 #endif
 399
 400 #if GMX_HAVE_SIMD_UPDATE
 401
 402 /*! \brief Load (aligned) the contents of GMX_SIMD_REAL_WIDTH rvec elements sequentially into 3 SIMD registers
 403  *
 404  * The loaded output is:
 405  * \p r0: { r[index][XX], r[index][YY], ... }
 406  * \p r1: { ... }
 407  * \p r2: { ..., r[index+GMX_SIMD_REAL_WIDTH-1][YY], r[index+GMX_SIMD_REAL_WIDTH-1][ZZ] }
 408  *
 409  * \param[in]  r      Real 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 load
 411  * \param[out] r0     Pointer to first SIMD register
 412  * \param[out] r1     Pointer to second SIMD register
 413  * \param[out] r2     Pointer to third SIMD register
 414  */
 415 static inline void simdLoadRvecs(const rvec* r, int index, SimdReal* r0, SimdReal* r1, SimdReal* r2)
 416 {
 417     const real* realPtr = r[index];
 418
 419     GMX_ASSERT(isSimdAligned(realPtr), "Pointer should be SIMD aligned");
 420
 421     *r0 = simdLoad(realPtr + 0 * GMX_SIMD_REAL_WIDTH);
 422     *r1 = simdLoad(realPtr + 1 * GMX_SIMD_REAL_WIDTH);
 423     *r2 = simdLoad(realPtr + 2 * GMX_SIMD_REAL_WIDTH);
 424 }
 425
 426 /*! \brief Store (aligned) 3 SIMD registers sequentially to GMX_SIMD_REAL_WIDTH rvec elements
 427  *
 428  * The stored output is:
 429  * \p r[index] = { { r0[0], r0[1], ... }
 430  * ...
 431  * \p r[index+GMX_SIMD_REAL_WIDTH-1] =  { ... , r2[GMX_SIMD_REAL_WIDTH-2], r2[GMX_SIMD_REAL_WIDTH-1] }
 432  *
 433  * \param[out] r      Pointer to an rvec array, has to be aligned to SIMD register width
 434  * \param[in]  index  Index of the first rvec triplet of reals to store to
 435  * \param[in]  r0     First SIMD register
 436  * \param[in]  r1     Second SIMD register
 437  * \param[in]  r2     Third SIMD register
 438  */
 439 static inline void simdStoreRvecs(rvec* r, int index, SimdReal r0, SimdReal r1, SimdReal r2)
 440 {
 441     real* realPtr = r[index];
 442
 443     GMX_ASSERT(isSimdAligned(realPtr), "Pointer should be SIMD aligned");
 444
 445     store(realPtr + 0 * GMX_SIMD_REAL_WIDTH, r0);
 446     store(realPtr + 1 * GMX_SIMD_REAL_WIDTH, r1);
 447     store(realPtr + 2 * GMX_SIMD_REAL_WIDTH, r2);
 448 }
 449
 450 /*! \brief Integrate using leap-frog with single group T-scaling and SIMD
 451  *
 452  * \tparam       storeUpdatedVelocities Tells whether we should store the updated velocities
 453  * \param[in]    start                  Index of first atom to update
 454  * \param[in]    nrend                  Last atom to update: \p nrend - 1
 455  * \param[in]    dt                     The time step
 456  * \param[in]    invMass                1/mass per atom
 457  * \param[in]    tcstat                 Temperature coupling information
 458  * \param[in]    x                      Input coordinates
 459  * \param[out]   xprime                 Updated coordinates
 460  * \param[inout] v                      Velocities, type either rvec* or const rvec*
 461  * \param[in]    f                      Forces
 462  */
 463 template<StoreUpdatedVelocities storeUpdatedVelocities, typename VelocityType>
 464 static std::enable_if_t<std::is_same<VelocityType, rvec>::value || std::is_same<VelocityType, const rvec>::value, void>
 465 updateMDLeapfrogSimpleSimd(int                               start,
 466                            int                               nrend,
 467                            real                              dt,
 468                            gmx::ArrayRef<const real>         invMass,
 469                            gmx::ArrayRef<const t_grp_tcstat> tcstat,
 470                            const rvec* gmx_restrict          x,
 471                            rvec* gmx_restrict                xprime,
 472                            VelocityType* gmx_restrict        v,
 473                            const rvec* gmx_restrict          f)
 474 {
 475     SimdReal timestep(dt);
 476     SimdReal lambdaSystem(tcstat[0].lambda);
 477
 478     /* We declare variables here, since code is often slower when declaring them inside the loop */
 479
 480     /* Note: We should implement a proper PaddedVector, so we don't need this check */
 481     GMX_ASSERT(isSimdAligned(invMass.data()), "invMass should be aligned");
 482
 483     for (int a = start; a < nrend; a += GMX_SIMD_REAL_WIDTH)
 484     {
 485         SimdReal invMass0, invMass1, invMass2;
 486         expandScalarsToTriplets(simdLoad(invMass.data() + a), &invMass0, &invMass1, &invMass2);
 487
 488         SimdReal v0, v1, v2;
 489         SimdReal f0, f1, f2;
 490         simdLoadRvecs(v, a, &v0, &v1, &v2);
 491         simdLoadRvecs(f, a, &f0, &f1, &f2);
 492
 493         v0 = fma(f0 * invMass0, timestep, lambdaSystem * v0);
 494         v1 = fma(f1 * invMass1, timestep, lambdaSystem * v1);
 495         v2 = fma(f2 * invMass2, timestep, lambdaSystem * v2);
 496
 497         if constexpr (storeUpdatedVelocities == StoreUpdatedVelocities::yes)
 498         {
 499             simdStoreRvecs(v, a, v0, v1, v2);
 500         }
 501
 502         // NOLINTNEXTLINE(readability-misleading-indentation) remove when clang-tidy-13 is required
 503         SimdReal x0, x1, x2;
 504         simdLoadRvecs(x, a, &x0, &x1, &x2);
 505
 506         SimdReal xprime0 = fma(v0, timestep, x0);
 507         SimdReal xprime1 = fma(v1, timestep, x1);
 508         SimdReal xprime2 = fma(v2, timestep, x2);
 509
 510         simdStoreRvecs(xprime, a, xprime0, xprime1, xprime2);
 511     }
 512 }
 513
 514 #endif // GMX_HAVE_SIMD_UPDATE
 515
 516 /*! \brief Sets the NEMD acceleration type */
 517 enum class AccelerationType
 518 {
 519     none,
 520     cosine
 521 };
 522
 523 /*! \brief Integrate using leap-frog with support for everything.
 524  *
 525  * \tparam        accelerationType  Type of NEMD acceleration.
 526  * \param[in]     start             Index of first atom to update.
 527  * \param[in]     nrend             Last atom to update: \p nrend - 1.
 528  * \param[in]     doNoseHoover      If to apply Nose-Hoover T-coupling.
 529  * \param[in]     dt                The time step.
 530  * \param[in]     dtPressureCouple  Time step for pressure coupling, is 0 when no pressure
 531  *                                  coupling should be applied at this step.
 532  * \param[in]     cTC               Temperature coupling group indices
 533  * \param[in]     invMassPerDim     Inverse mass per dimension
 534  * \param[in]     ekind             Kinetic energy data.
 535  * \param[in]     box               The box dimensions.
 536  * \param[in]     x                 Input coordinates.
 537  * \param[out]    xprime            Updated coordinates.
 538  * \param[inout]  v                 Velocities.
 539  * \param[in]     f                 Forces.
 540  * \param[in]     nh_vxi            Nose-Hoover velocity scaling factors.
 541  * \param[in]     nsttcouple        Frequency of the temperature coupling steps.
 542  * \param[in]     M                 Parrinello-Rahman scaling matrix.
 543  */
 544 template<AccelerationType accelerationType>
 545 static void updateMDLeapfrogGeneral(int                                 start,
 546                                     int                                 nrend,
 547                                     bool                                doNoseHoover,
 548                                     real                                dt,
 549                                     real                                dtPressureCouple,
 550                                     gmx::ArrayRef<const unsigned short> cTC,
 551                                     gmx::ArrayRef<const rvec>           invMassPerDim,
 552                                     const gmx_ekindata_t*               ekind,
 553                                     const matrix                        box,
 554                                     const rvec* gmx_restrict            x,
 555                                     rvec* gmx_restrict                  xprime,
 556                                     rvec* gmx_restrict                  v,
 557                                     const rvec* gmx_restrict            f,
 558                                     const double* gmx_restrict          nh_vxi,
 559                                     const int                           nsttcouple,
 560                                     const matrix                        M)
 561 {
 562     /* This is a version of the leap-frog integrator that supports
 563      * all combinations of T-coupling, P-coupling and NEMD.
 564      * Nose-Hoover uses the reversible leap-frog integrator from
 565      * Holian et al. Phys Rev E 52(3) : 2338, 1995
 566      */
 567
 568     gmx::ArrayRef<const t_grp_tcstat> tcstat = ekind->tcstat;
 569
 570     /* Initialize group values, changed later when multiple groups are used */
 571     int gt = 0;
 572
 573     real omega_Z = 2 * static_cast<real>(M_PI) / box[ZZ][ZZ];
 574
 575     for (int n = start; n < nrend; n++)
 576     {
 577         if (!cTC.empty())
 578         {
 579             gt = cTC[n];
 580         }
 581         real lg = tcstat[gt].lambda;
 582
 583         rvec vRel;
 584         real cosineZ, vCosine;
 585         switch (accelerationType)
 586         {
 587             case AccelerationType::none: copy_rvec(v[n], vRel); break;
 588             case AccelerationType::cosine:
 589                 cosineZ = std::cos(x[n][ZZ] * omega_Z);
 590                 vCosine = cosineZ * ekind->cosacc.vcos;
 591                 /* Avoid scaling the cosine profile velocity */
 592                 copy_rvec(v[n], vRel);
 593                 vRel[XX] -= vCosine;
 594                 break;
 595         }
 596
 597         real factorNH = 0.0;
 598         if (doNoseHoover)
 599         {
 600             /* Here we account for multiple time stepping, by increasing
 601              * the Nose-Hoover correction by nsttcouple
 602              * TODO: This can be pre-computed.
 603              */
 604             factorNH = 0.5 * nsttcouple * dt * nh_vxi[gt];
 605         }
 606
 607         for (int d = 0; d < DIM; d++)
 608         {
 609             real vNew = (lg * vRel[d]
 610                          + (f[n][d] * invMassPerDim[n][d] * dt - factorNH * vRel[d]
 611                             - dtPressureCouple * iprod(M[d], vRel)))
 612                         / (1 + factorNH);
 613             switch (accelerationType)
 614             {
 615                 case AccelerationType::none: break;
 616                 case AccelerationType::cosine:
 617                     if (d == XX)
 618                     {
 619                         /* Add back the mean velocity and apply acceleration */
 620                         vNew += vCosine + cosineZ * ekind->cosacc.cos_accel * dt;
 621                     }
 622                     break;
 623             }
 624             v[n][d]      = vNew;
 625             xprime[n][d] = x[n][d] + vNew * dt;
 626         }
 627     }
 628 }
 629
 630 /*! \brief Handles the Leap-frog MD x and v integration */
 631 static void do_update_md(int                                  start,
 632                          int                                  nrend,
 633                          real                                 dt,
 634                          int64_t                              step,
 635                          const rvec* gmx_restrict             x,
 636                          rvec* gmx_restrict                   xprime,
 637                          rvec* gmx_restrict                   v,
 638                          const rvec* gmx_restrict             f,
 639                          const TemperatureCoupling            etc,
 640                          const PressureCoupling               epc,
 641                          const int                            nsttcouple,
 642                          const int                            nstpcouple,
 643                          gmx::ArrayRef<const unsigned short>  cTC,
 644                          gmx::ArrayRef<const real> gmx_unused invmass,
 645                          gmx::ArrayRef<const rvec>            invMassPerDim,
 646                          const gmx_ekindata_t*                ekind,
 647                          const matrix                         box,
 648                          const double* gmx_restrict           nh_vxi,
 649                          const matrix                         M,
 650                          bool gmx_unused                      havePartiallyFrozenAtoms)
 651 {
 652     GMX_ASSERT(nrend == start || xprime != x,
 653                "For SIMD optimization certain compilers need to have xprime != x");
 654
 655     /* Note: Berendsen pressure scaling is handled after do_update_md() */
 656     bool doTempCouple =
 657             (etc != TemperatureCoupling::No && do_per_step(step + nsttcouple - 1, nsttcouple));
 658     bool doNoseHoover       = (etc == TemperatureCoupling::NoseHoover && doTempCouple);
 659     bool doParrinelloRahman = (epc == PressureCoupling::ParrinelloRahman
 660                                && do_per_step(step + nstpcouple - 1, nstpcouple));
 661     bool doPROffDiagonal = (doParrinelloRahman && (M[YY][XX] != 0 || M[ZZ][XX] != 0 || M[ZZ][YY] != 0));
 662
 663     real dtPressureCouple = (doParrinelloRahman ? nstpcouple * dt : 0);
 664
 665     /* NEMD cosine acceleration is applied in updateMDLeapFrogGeneral */
 666     bool doAcceleration = (ekind->cosacc.cos_accel != 0);
 667
 668     if (doNoseHoover || doPROffDiagonal || doAcceleration)
 669     {
 670         matrix stepM;
 671         if (!doParrinelloRahman)
 672         {
 673             /* We should not apply PR scaling at this step */
 674             clear_mat(stepM);
 675         }
 676         else
 677         {
 678             copy_mat(M, stepM);
 679         }
 680
 681         if (!doAcceleration)
 682         {
 683             updateMDLeapfrogGeneral<AccelerationType::none>(start,
 684                                                             nrend,
 685                                                             doNoseHoover,
 686                                                             dt,
 687                                                             dtPressureCouple,
 688                                                             cTC,
 689                                                             invMassPerDim,
 690                                                             ekind,
 691                                                             box,
 692                                                             x,
 693                                                             xprime,
 694                                                             v,
 695                                                             f,
 696                                                             nh_vxi,
 697                                                             nsttcouple,
 698                                                             stepM);
 699         }
 700         else
 701         {
 702             updateMDLeapfrogGeneral<AccelerationType::cosine>(start,
 703                                                               nrend,
 704                                                               doNoseHoover,
 705                                                               dt,
 706                                                               dtPressureCouple,
 707                                                               cTC,
 708                                                               invMassPerDim,
 709                                                               ekind,
 710                                                               box,
 711                                                               x,
 712                                                               xprime,
 713                                                               v,
 714                                                               f,
 715                                                               nh_vxi,
 716                                                               nsttcouple,
 717                                                               stepM);
 718         }
 719     }
 720     else
 721     {
 722         /* We determine if we have a single T-coupling lambda value for all
 723          * atoms. That allows for better SIMD acceleration in the template.
 724          * If we do not do temperature coupling (in the run or this step),
 725          * all scaling values are 1, so we effectively have a single value.
 726          */
 727         bool haveSingleTempScaleValue = (!doTempCouple || ekind->ngtc == 1);
 728
 729         /* Extract some pointers needed by all cases */
 730         gmx::ArrayRef<const t_grp_tcstat> tcstat = ekind->tcstat;
 731
 732         if (doParrinelloRahman)
 733         {
 734             GMX_ASSERT(!doPROffDiagonal,
 735                        "updateMDLeapfrogSimple only support diagonal Parrinello-Rahman scaling "
 736                        "matrices");
 737
 738             rvec diagM;
 739             for (int d = 0; d < DIM; d++)
 740             {
 741                 diagM[d] = M[d][d];
 742             }
 743
 744             if (haveSingleTempScaleValue)
 745             {
 746                 updateMDLeapfrogSimple<StoreUpdatedVelocities::yes, NumTempScaleValues::single, ApplyParrinelloRahmanVScaling::diagonal>(
 747                         start, nrend, dt, dtPressureCouple, invMassPerDim, tcstat, cTC, diagM, x, xprime, v, f);
 748             }
 749             else
 750             {
 751                 updateMDLeapfrogSimple<StoreUpdatedVelocities::yes, NumTempScaleValues::multiple, ApplyParrinelloRahmanVScaling::diagonal>(
 752                         start, nrend, dt, dtPressureCouple, invMassPerDim, tcstat, cTC, diagM, x, xprime, v, f);
 753             }
 754         }
 755         else
 756         {
 757             if (haveSingleTempScaleValue)
 758             {
 759                 /* Note that modern compilers are pretty good at vectorizing
 760                  * updateMDLeapfrogSimple(). But the SIMD version will still
 761                  * be faster because invMass lowers the cache pressure
 762                  * compared to invMassPerDim.
 763                  */
 764 #if GMX_HAVE_SIMD_UPDATE
 765                 /* Check if we can use invmass instead of invMassPerDim */
 766                 if (!havePartiallyFrozenAtoms)
 767                 {
 768                     updateMDLeapfrogSimpleSimd<StoreUpdatedVelocities::yes>(
 769                             start, nrend, dt, invmass, tcstat, x, xprime, v, f);
 770                 }
 771                 else
 772 #endif
 773                 {
 774                     updateMDLeapfrogSimple<StoreUpdatedVelocities::yes, NumTempScaleValues::single, ApplyParrinelloRahmanVScaling::no>(
 775                             start, nrend, dt, dtPressureCouple, invMassPerDim, tcstat, cTC, nullptr, x, xprime, v, f);
 776                 }
 777             }
 778             else
 779             {
 780                 updateMDLeapfrogSimple<StoreUpdatedVelocities::yes, NumTempScaleValues::multiple, ApplyParrinelloRahmanVScaling::no>(
 781                         start, nrend, dt, dtPressureCouple, invMassPerDim, tcstat, cTC, nullptr, x, xprime, v, f);
 782             }
 783         }
 784     }
 785 }
 786 /*! \brief Handles the Leap-frog MD x and v integration */
 787 static void doUpdateMDDoNotUpdateVelocities(int                      start,
 788                                             int                      nrend,
 789                                             real                     dt,
 790                                             const rvec* gmx_restrict x,
 791                                             rvec* gmx_restrict       xprime,
 792                                             const rvec* gmx_restrict v,
 793                                             const rvec* gmx_restrict f,
 794                                             bool gmx_unused          havePartiallyFrozenAtoms,
 795                                             gmx::ArrayRef<const real> gmx_unused invmass,
 796                                             gmx::ArrayRef<const rvec>            invMassPerDim,
 797                                             const gmx_ekindata_t&                ekind)
 798 {
 799     GMX_ASSERT(nrend == start || xprime != x,
 800                "For SIMD optimization certain compilers need to have xprime != x");
 801
 802     gmx::ArrayRef<const t_grp_tcstat> tcstat = ekind.tcstat;
 803
 804     /* Check if we can use invmass instead of invMassPerDim */
 805 #if GMX_HAVE_SIMD_UPDATE
 806     if (!havePartiallyFrozenAtoms)
 807     {
 808         updateMDLeapfrogSimpleSimd<StoreUpdatedVelocities::no>(
 809                 start, nrend, dt, invmass, tcstat, x, xprime, v, f);
 810     }
 811     else
 812 #endif
 813     {
 814         updateMDLeapfrogSimple<StoreUpdatedVelocities::no, NumTempScaleValues::single, ApplyParrinelloRahmanVScaling::no>(
 815                 start, nrend, dt, dt, invMassPerDim, tcstat, gmx::ArrayRef<const unsigned short>(), nullptr, x, xprime, v, f);
 816     }
 817 }
 818
 819 static void do_update_vv_vel(int                                 start,
 820                              int                                 nrend,
 821                              real                                dt,
 822                              gmx::ArrayRef<const ivec>           nFreeze,
 823                              gmx::ArrayRef<const real>           invmass,
 824                              gmx::ArrayRef<const ParticleType>   ptype,
 825                              gmx::ArrayRef<const unsigned short> cFREEZE,
 826                              rvec                                v[],
 827                              const rvec                          f[],
 828                              gmx_bool                            bExtended,
 829                              real                                veta,
 830                              real                                alpha)
 831 {
 832     int  gf = 0;
 833     int  n, d;
 834     real g, mv1, mv2;
 835
 836     if (bExtended)
 837     {
 838         g   = 0.25 * dt * veta * alpha;
 839         mv1 = std::exp(-g);
 840         mv2 = gmx::series_sinhx(g);
 841     }
 842     else
 843     {
 844         mv1 = 1.0;
 845         mv2 = 1.0;
 846     }
 847     for (n = start; n < nrend; n++)
 848     {
 849         real w_dt = invmass[n] * dt;
 850         if (!cFREEZE.empty())
 851         {
 852             gf = cFREEZE[n];
 853         }
 854
 855         for (d = 0; d < DIM; d++)
 856         {
 857             if ((ptype[n] != ParticleType::Shell) && !nFreeze[gf][d])
 858             {
 859                 v[n][d] = mv1 * (mv1 * v[n][d] + 0.5 * (w_dt * mv2 * f[n][d]));
 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                              gmx::ArrayRef<const ivec>           nFreeze,
 873                              gmx::ArrayRef<const ParticleType>   ptype,
 874                              gmx::ArrayRef<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.empty())
 902         {
 903             gf = cFREEZE[n];
 904         }
 905
 906         for (d = 0; d < DIM; d++)
 907         {
 908             if ((ptype[n] != ParticleType::Shell) && !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 == IntegrationAlgorithm::BD)
 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]        = gmx::c_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 == IntegrationAlgorithm::BD)
 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 * gmx::c_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 * gmx::c_boltz * inputRecord.opts.ref_t[gt]);
 989             }
 990         }
 991     }
 992     if (inputRecord.eI == IntegrationAlgorithm::SD1)
 993     {
 994         for (int gt = 0; gt < inputRecord.opts.ngtc; gt++)
 995         {
 996             real kT = gmx::c_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), deform_(boxDeformation)
1005 {
1006     update_temperature_constants(inputRecord);
1007     xp_.resizeWithPadding(0);
1008 }
1009
1010 void Update::updateAfterPartition(int                                 numAtoms,
1011                                   gmx::ArrayRef<const unsigned short> cFREEZE,
1012                                   gmx::ArrayRef<const unsigned short> cTC)
1013 {
1014
1015     impl_->xp()->resizeWithPadding(numAtoms);
1016     impl_->cFREEZE_ = cFREEZE;
1017     impl_->cTC_     = cTC;
1018 }
1019
1020 /*! \brief Sets the SD update type */
1021 enum class SDUpdate : int
1022 {
1023     ForcesOnly,
1024     FrictionAndNoiseOnly,
1025     Combined
1026 };
1027
1028 /*! \brief SD integrator update
1029  *
1030  * Two phases are required in the general case of a constrained
1031  * update, the first phase from the contribution of forces, before
1032  * applying constraints, and then a second phase applying the friction
1033  * and noise, and then further constraining. For details, see
1034  * Goga2012.
1035  *
1036  * Without constraints, the two phases can be combined, for
1037  * efficiency.
1038  *
1039  * Thus three instantiations of this templated function will be made,
1040  * two with only one contribution, and one with both contributions. */
1041 template<SDUpdate updateType>
1042 static void doSDUpdateGeneral(const gmx_stochd_t&                 sd,
1043                               int                                 start,
1044                               int                                 nrend,
1045                               real                                dt,
1046                               gmx::ArrayRef<const ivec>           nFreeze,
1047                               gmx::ArrayRef<const real>           invmass,
1048                               gmx::ArrayRef<const ParticleType>   ptype,
1049                               gmx::ArrayRef<const unsigned short> cFREEZE,
1050                               gmx::ArrayRef<const unsigned short> cTC,
1051                               const rvec                          x[],
1052                               rvec                                xprime[],
1053                               rvec                                v[],
1054                               const rvec                          f[],
1055                               int64_t                             step,
1056                               int                                 seed,
1057                               const int*                          gatindex)
1058 {
1059     // cTC and cFREEZE can be nullptr any time, but various
1060     // instantiations do not make sense with particular pointer
1061     // values.
1062     if (updateType == SDUpdate::ForcesOnly)
1063     {
1064         GMX_ASSERT(f != nullptr, "SD update with only forces requires forces");
1065         GMX_ASSERT(cTC.empty(), "SD update with only forces cannot handle temperature groups");
1066     }
1067     if (updateType == SDUpdate::FrictionAndNoiseOnly)
1068     {
1069         GMX_ASSERT(f == nullptr, "SD update with only noise cannot handle forces");
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.empty() ? cFREEZE[n] : 0;
1090         int temperatureGroup = !cTC.empty() ? cTC[n] : 0;
1091
1092         for (int d = 0; d < DIM; d++)
1093         {
1094             if ((ptype[n] != ParticleType::Shell) && !nFreeze[freezeGroup][d])
1095             {
1096                 if (updateType == SDUpdate::ForcesOnly)
1097                 {
1098                     real vn = v[n][d] + inverseMass * f[n][d] * dt;
1099                     v[n][d] = vn;
1100                     // Simple position update.
1101                     xprime[n][d] = x[n][d] + v[n][d] * dt;
1102                 }
1103                 else if (updateType == SDUpdate::FrictionAndNoiseOnly)
1104                 {
1105                     real vn = v[n][d];
1106                     v[n][d] = (vn * sd.sdc[temperatureGroup].em
1107                                + invsqrtMass * sd.sdsig[temperatureGroup].V * dist(rng));
1108                     // The previous phase already updated the
1109                     // positions with a full v*dt term that must
1110                     // now be half removed.
1111                     xprime[n][d] = xprime[n][d] + 0.5 * (v[n][d] - vn) * dt;
1112                 }
1113                 else
1114                 {
1115                     real vn = v[n][d] + inverseMass * f[n][d] * dt;
1116                     v[n][d] = (vn * sd.sdc[temperatureGroup].em
1117                                + invsqrtMass * sd.sdsig[temperatureGroup].V * dist(rng));
1118                     // Here we include half of the friction+noise
1119                     // update of v into the position update.
1120                     xprime[n][d] = x[n][d] + 0.5 * (vn + v[n][d]) * dt;
1121                 }
1122             }
1123             else
1124             {
1125                 // When using constraints, the update is split into
1126                 // two phases, but we only need to zero the update of
1127                 // virtual, shell or frozen particles in at most one
1128                 // of the phases.
1129                 if (updateType != SDUpdate::FrictionAndNoiseOnly)
1130                 {
1131                     v[n][d]      = 0.0;
1132                     xprime[n][d] = x[n][d];
1133                 }
1134             }
1135         }
1136     }
1137 }
1138
1139 static void do_update_sd(int                                 start,
1140                          int                                 nrend,
1141                          real                                dt,
1142                          int64_t                             step,
1143                          const rvec* gmx_restrict            x,
1144                          rvec* gmx_restrict                  xprime,
1145                          rvec* gmx_restrict                  v,
1146                          const rvec* gmx_restrict            f,
1147                          gmx::ArrayRef<const ivec>           nFreeze,
1148                          gmx::ArrayRef<const real>           invmass,
1149                          gmx::ArrayRef<const ParticleType>   ptype,
1150                          gmx::ArrayRef<const unsigned short> cFREEZE,
1151                          gmx::ArrayRef<const unsigned short> cTC,
1152                          int                                 seed,
1153                          const t_commrec*                    cr,
1154                          const gmx_stochd_t&                 sd,
1155                          bool                                haveConstraints)
1156 {
1157     if (haveConstraints)
1158     {
1159         // With constraints, the SD update is done in 2 parts
1160         doSDUpdateGeneral<SDUpdate::ForcesOnly>(sd,
1161                                                 start,
1162                                                 nrend,
1163                                                 dt,
1164                                                 nFreeze,
1165                                                 invmass,
1166                                                 ptype,
1167                                                 cFREEZE,
1168                                                 gmx::ArrayRef<const unsigned short>(),
1169                                                 x,
1170                                                 xprime,
1171                                                 v,
1172                                                 f,
1173                                                 step,
1174                                                 seed,
1175                                                 nullptr);
1176     }
1177     else
1178     {
1179         doSDUpdateGeneral<SDUpdate::Combined>(sd,
1180                                               start,
1181                                               nrend,
1182                                               dt,
1183                                               nFreeze,
1184                                               invmass,
1185                                               ptype,
1186                                               cFREEZE,
1187                                               cTC,
1188                                               x,
1189                                               xprime,
1190                                               v,
1191                                               f,
1192                                               step,
1193                                               seed,
1194                                               DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1195     }
1196 }
1197
1198 static void do_update_bd(int                                 start,
1199                          int                                 nrend,
1200                          real                                dt,
1201                          int64_t                             step,
1202                          const rvec* gmx_restrict            x,
1203                          rvec* gmx_restrict                  xprime,
1204                          rvec* gmx_restrict                  v,
1205                          const rvec* gmx_restrict            f,
1206                          gmx::ArrayRef<const ivec>           nFreeze,
1207                          gmx::ArrayRef<const real>           invmass,
1208                          gmx::ArrayRef<const ParticleType>   ptype,
1209                          gmx::ArrayRef<const unsigned short> cFREEZE,
1210                          gmx::ArrayRef<const unsigned short> cTC,
1211                          real                                friction_coefficient,
1212                          const real*                         rf,
1213                          int                                 seed,
1214                          const int*                          gatindex)
1215 {
1216     /* note -- these appear to be full step velocities . . .  */
1217     int  gf = 0, gt = 0;
1218     real vn;
1219     real invfr = 0;
1220     int  n, d;
1221     // Use 1 bit of internal counters to give us 2*2 64-bits values per stream
1222     // Each 64-bit value is enough for 4 normal distribution table numbers.
1223     gmx::ThreeFry2x64<0>                       rng(seed, gmx::RandomDomain::UpdateCoordinates);
1224     gmx::TabulatedNormalDistribution<real, 14> dist;
1225
1226     if (friction_coefficient != 0)
1227     {
1228         invfr = 1.0 / friction_coefficient;
1229     }
1230
1231     for (n = start; (n < nrend); n++)
1232     {
1233         int ng = gatindex ? gatindex[n] : n;
1234
1235         rng.restart(step, ng);
1236         dist.reset();
1237
1238         if (!cFREEZE.empty())
1239         {
1240             gf = cFREEZE[n];
1241         }
1242         if (!cTC.empty())
1243         {
1244             gt = cTC[n];
1245         }
1246         for (d = 0; (d < DIM); d++)
1247         {
1248             if ((ptype[n] != ParticleType::Shell) && !nFreeze[gf][d])
1249             {
1250                 if (friction_coefficient != 0)
1251                 {
1252                     vn = invfr * f[n][d] + rf[gt] * dist(rng);
1253                 }
1254                 else
1255                 {
1256                     /* NOTE: invmass = 2/(mass*friction_constant*dt) */
1257                     vn = 0.5 * invmass[n] * f[n][d] * dt
1258                          + std::sqrt(0.5 * invmass[n]) * rf[gt] * dist(rng);
1259                 }
1260
1261                 v[n][d]      = vn;
1262                 xprime[n][d] = x[n][d] + vn * dt;
1263             }
1264             else
1265             {
1266                 v[n][d]      = 0.0;
1267                 xprime[n][d] = x[n][d];
1268             }
1269         }
1270     }
1271 }
1272
1273 extern void init_ekinstate(ekinstate_t* ekinstate, const t_inputrec* ir)
1274 {
1275     ekinstate->ekin_n = ir->opts.ngtc;
1276     snew(ekinstate->ekinh, ekinstate->ekin_n);
1277     snew(ekinstate->ekinf, ekinstate->ekin_n);
1278     snew(ekinstate->ekinh_old, ekinstate->ekin_n);
1279     ekinstate->ekinscalef_nhc.resize(ekinstate->ekin_n);
1280     ekinstate->ekinscaleh_nhc.resize(ekinstate->ekin_n);
1281     ekinstate->vscale_nhc.resize(ekinstate->ekin_n);
1282     ekinstate->dekindl          = 0;
1283     ekinstate->mvcos            = 0;
1284     ekinstate->hasReadEkinState = false;
1285 }
1286
1287 void update_ekinstate(ekinstate_t* ekinstate, const gmx_ekindata_t* ekind)
1288 {
1289     int i;
1290
1291     for (i = 0; i < ekinstate->ekin_n; i++)
1292     {
1293         copy_mat(ekind->tcstat[i].ekinh, ekinstate->ekinh[i]);
1294         copy_mat(ekind->tcstat[i].ekinf, ekinstate->ekinf[i]);
1295         copy_mat(ekind->tcstat[i].ekinh_old, ekinstate->ekinh_old[i]);
1296         ekinstate->ekinscalef_nhc[i] = ekind->tcstat[i].ekinscalef_nhc;
1297         ekinstate->ekinscaleh_nhc[i] = ekind->tcstat[i].ekinscaleh_nhc;
1298         ekinstate->vscale_nhc[i]     = ekind->tcstat[i].vscale_nhc;
1299     }
1300
1301     copy_mat(ekind->ekin, ekinstate->ekin_total);
1302     ekinstate->dekindl = ekind->dekindl;
1303     ekinstate->mvcos   = ekind->cosacc.mvcos;
1304 }
1305
1306 void restore_ekinstate_from_state(const t_commrec* cr, gmx_ekindata_t* ekind, const ekinstate_t* ekinstate)
1307 {
1308     int i, n;
1309
1310     if (MASTER(cr))
1311     {
1312         for (i = 0; i < ekinstate->ekin_n; i++)
1313         {
1314             copy_mat(ekinstate->ekinh[i], ekind->tcstat[i].ekinh);
1315             copy_mat(ekinstate->ekinf[i], ekind->tcstat[i].ekinf);
1316             copy_mat(ekinstate->ekinh_old[i], ekind->tcstat[i].ekinh_old);
1317             ekind->tcstat[i].ekinscalef_nhc = ekinstate->ekinscalef_nhc[i];
1318             ekind->tcstat[i].ekinscaleh_nhc = ekinstate->ekinscaleh_nhc[i];
1319             ekind->tcstat[i].vscale_nhc     = ekinstate->vscale_nhc[i];
1320         }
1321
1322         copy_mat(ekinstate->ekin_total, ekind->ekin);
1323
1324         ekind->dekindl      = ekinstate->dekindl;
1325         ekind->cosacc.mvcos = ekinstate->mvcos;
1326         n                   = ekinstate->ekin_n;
1327     }
1328
1329     if (PAR(cr))
1330     {
1331         gmx_bcast(sizeof(n), &n, cr->mpi_comm_mygroup);
1332         for (i = 0; i < n; i++)
1333         {
1334             gmx_bcast(DIM * DIM * sizeof(ekind->tcstat[i].ekinh[0][0]),
1335                       ekind->tcstat[i].ekinh[0],
1336                       cr->mpi_comm_mygroup);
1337             gmx_bcast(DIM * DIM * sizeof(ekind->tcstat[i].ekinf[0][0]),
1338                       ekind->tcstat[i].ekinf[0],
1339                       cr->mpi_comm_mygroup);
1340             gmx_bcast(DIM * DIM * sizeof(ekind->tcstat[i].ekinh_old[0][0]),
1341                       ekind->tcstat[i].ekinh_old[0],
1342                       cr->mpi_comm_mygroup);
1343
1344             gmx_bcast(sizeof(ekind->tcstat[i].ekinscalef_nhc),
1345                       &(ekind->tcstat[i].ekinscalef_nhc),
1346                       cr->mpi_comm_mygroup);
1347             gmx_bcast(sizeof(ekind->tcstat[i].ekinscaleh_nhc),
1348                       &(ekind->tcstat[i].ekinscaleh_nhc),
1349                       cr->mpi_comm_mygroup);
1350             gmx_bcast(sizeof(ekind->tcstat[i].vscale_nhc), &(ekind->tcstat[i].vscale_nhc), cr->mpi_comm_mygroup);
1351         }
1352         gmx_bcast(DIM * DIM * sizeof(ekind->ekin[0][0]), ekind->ekin[0], cr->mpi_comm_mygroup);
1353
1354         gmx_bcast(sizeof(ekind->dekindl), &ekind->dekindl, cr->mpi_comm_mygroup);
1355         gmx_bcast(sizeof(ekind->cosacc.mvcos), &ekind->cosacc.mvcos, cr->mpi_comm_mygroup);
1356     }
1357 }
1358
1359 void getThreadAtomRange(int numThreads, int threadIndex, int numAtoms, int* startAtom, int* endAtom)
1360 {
1361 #if GMX_HAVE_SIMD_UPDATE
1362     constexpr int blockSize = GMX_SIMD_REAL_WIDTH;
1363 #else
1364     constexpr int blockSize = 1;
1365 #endif
1366     int numBlocks = (numAtoms + blockSize - 1) / blockSize;
1367
1368     *startAtom = ((numBlocks * threadIndex) / numThreads) * blockSize;
1369     *endAtom   = ((numBlocks * (threadIndex + 1)) / numThreads) * blockSize;
1370     if (threadIndex == numThreads - 1)
1371     {
1372         *endAtom = numAtoms;
1373     }
1374 }
1375
1376 void Update::Impl::update_sd_second_half(const t_inputrec&                 inputRecord,
1377                                          int64_t                           step,
1378                                          real*                             dvdlambda,
1379                                          int                               homenr,
1380                                          gmx::ArrayRef<const ParticleType> ptype,
1381                                          gmx::ArrayRef<const real>         invMass,
1382                                          t_state*                          state,
1383                                          const t_commrec*                  cr,
1384                                          t_nrnb*                           nrnb,
1385                                          gmx_wallcycle*                    wcycle,
1386                                          gmx::Constraints*                 constr,
1387                                          bool                              do_log,
1388                                          bool                              do_ene)
1389 {
1390     if (!constr)
1391     {
1392         return;
1393     }
1394     if (inputRecord.eI == IntegrationAlgorithm::SD1)
1395     {
1396
1397         /* Cast delta_t from double to real to make the integrators faster.
1398          * The only reason for having delta_t double is to get accurate values
1399          * for t=delta_t*step when step is larger than float precision.
1400          * For integration dt the accuracy of real suffices, since with
1401          * integral += dt*integrand the increment is nearly always (much) smaller
1402          * than the integral (and the integrand has real precision).
1403          */
1404         real dt = inputRecord.delta_t;
1405
1406         wallcycle_start(wcycle, WallCycleCounter::Update);
1407
1408         int nth = gmx_omp_nthreads_get(ModuleMultiThread::Update);
1409
1410 #pragma omp parallel for num_threads(nth) schedule(static)
1411         for (int th = 0; th < nth; th++)
1412         {
1413             try
1414             {
1415                 int start_th, end_th;
1416                 getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
1417
1418                 doSDUpdateGeneral<SDUpdate::FrictionAndNoiseOnly>(
1419                         sd_,
1420                         start_th,
1421                         end_th,
1422                         dt,
1423                         gmx::arrayRefFromArray(inputRecord.opts.nFreeze, inputRecord.opts.ngfrz),
1424                         invMass,
1425                         ptype,
1426                         cFREEZE_,
1427                         cTC_,
1428                         state->x.rvec_array(),
1429                         xp_.rvec_array(),
1430                         state->v.rvec_array(),
1431                         nullptr,
1432                         step,
1433                         inputRecord.ld_seed,
1434                         DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1435             }
1436             GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1437         }
1438         inc_nrnb(nrnb, eNR_UPDATE, homenr);
1439         wallcycle_stop(wcycle, WallCycleCounter::Update);
1440
1441         /* Constrain the coordinates upd->xp for half a time step */
1442         bool computeVirial = false;
1443         constr->apply(do_log,
1444                       do_ene,
1445                       step,
1446                       1,
1447                       0.5,
1448                       state->x.arrayRefWithPadding(),
1449                       xp_.arrayRefWithPadding(),
1450                       ArrayRef<RVec>(),
1451                       state->box,
1452                       state->lambda[static_cast<int>(FreeEnergyPerturbationCouplingType::Bonded)],
1453                       dvdlambda,
1454                       state->v.arrayRefWithPadding(),
1455                       computeVirial,
1456                       nullptr,
1457                       ConstraintVariable::Positions);
1458     }
1459 }
1460
1461 void Update::Impl::finish_update(const t_inputrec&                   inputRecord,
1462                                  const bool                          havePartiallyFrozenAtoms,
1463                                  const int                           homenr,
1464                                  gmx::ArrayRef<const unsigned short> cFREEZE,
1465                                  t_state*                            state,
1466                                  gmx_wallcycle*                      wcycle,
1467                                  const bool                          haveConstraints)
1468 {
1469     /* NOTE: Currently we always integrate to a temporary buffer and
1470      * then copy the results back here.
1471      */
1472
1473     wallcycle_start_nocount(wcycle, WallCycleCounter::Update);
1474
1475     auto xp = makeConstArrayRef(xp_).subArray(0, homenr);
1476     auto x  = makeArrayRef(state->x).subArray(0, homenr);
1477
1478     if (havePartiallyFrozenAtoms && haveConstraints)
1479     {
1480         /* We have atoms that are frozen along some, but not all dimensions,
1481          * then constraints will have moved them also along the frozen dimensions.
1482          * To freeze such degrees of freedom we do not copy them back here.
1483          */
1484         const ivec* nFreeze = inputRecord.opts.nFreeze;
1485
1486         for (int i = 0; i < homenr; i++)
1487         {
1488             const int g = cFREEZE[i];
1489
1490             for (int d = 0; d < DIM; d++)
1491             {
1492                 if (nFreeze[g][d] == 0)
1493                 {
1494                     x[i][d] = xp[i][d];
1495                 }
1496             }
1497         }
1498     }
1499     else
1500     {
1501         /* We have no frozen atoms or fully frozen atoms which have not
1502          * been moved by the update, so we can simply copy all coordinates.
1503          */
1504         int gmx_unused nth = gmx_omp_nthreads_get(ModuleMultiThread::Update);
1505 #pragma omp parallel for num_threads(nth) schedule(static)
1506         for (int i = 0; i < homenr; i++)
1507         {
1508             // Trivial statement, does not throw
1509             x[i] = xp[i];
1510         }
1511     }
1512
1513     wallcycle_stop(wcycle, WallCycleCounter::Update);
1514 }
1515
1516 void Update::Impl::update_coords(const t_inputrec&                 inputRecord,
1517                                  int64_t                           step,
1518                                  int                               homenr,
1519                                  bool                              havePartiallyFrozenAtoms,
1520                                  gmx::ArrayRef<const ParticleType> ptype,
1521                                  gmx::ArrayRef<const real>         invMass,
1522                                  gmx::ArrayRef<const rvec>         invMassPerDim,
1523                                  t_state*                          state,
1524                                  const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
1525                                  t_fcdata*                                        fcdata,
1526                                  const gmx_ekindata_t*                            ekind,
1527                                  const matrix                                     M,
1528                                  int                                              updatePart,
1529                                  const t_commrec*                                 cr,
1530                                  const bool                                       haveConstraints)
1531 {
1532     /* Running the velocity half does nothing except for velocity verlet */
1533     if ((updatePart == etrtVELOCITY1 || updatePart == etrtVELOCITY2) && !EI_VV(inputRecord.eI))
1534     {
1535         gmx_incons("update_coords called for velocity without VV integrator");
1536     }
1537
1538     /* Cast to real for faster code, no loss in precision (see comment above) */
1539     real dt = inputRecord.delta_t;
1540
1541     /* We need to update the NMR restraint history when time averaging is used */
1542     if (state->flags & enumValueToBitMask(StateEntry::DisreRm3Tav))
1543     {
1544         update_disres_history(*fcdata->disres, &state->hist);
1545     }
1546     if (state->flags & enumValueToBitMask(StateEntry::OrireDtav))
1547     {
1548         GMX_ASSERT(fcdata, "Need valid fcdata");
1549         fcdata->orires->updateHistory();
1550     }
1551
1552     /* ############# START The update of velocities and positions ######### */
1553     int nth = gmx_omp_nthreads_get(ModuleMultiThread::Update);
1554
1555 #pragma omp parallel for num_threads(nth) schedule(static)
1556     for (int th = 0; th < nth; th++)
1557     {
1558         try
1559         {
1560             int start_th, end_th;
1561             getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
1562
1563             const rvec* x_rvec  = state->x.rvec_array();
1564             rvec*       xp_rvec = xp_.rvec_array();
1565             rvec*       v_rvec  = state->v.rvec_array();
1566             const rvec* f_rvec  = as_rvec_array(f.unpaddedConstArrayRef().data());
1567
1568             switch (inputRecord.eI)
1569             {
1570                 case (IntegrationAlgorithm::MD):
1571                     do_update_md(start_th,
1572                                  end_th,
1573                                  dt,
1574                                  step,
1575                                  x_rvec,
1576                                  xp_rvec,
1577                                  v_rvec,
1578                                  f_rvec,
1579                                  inputRecord.etc,
1580                                  inputRecord.epc,
1581                                  inputRecord.nsttcouple,
1582                                  inputRecord.nstpcouple,
1583                                  cTC_,
1584                                  invMass,
1585                                  invMassPerDim,
1586                                  ekind,
1587                                  state->box,
1588                                  state->nosehoover_vxi.data(),
1589                                  M,
1590                                  havePartiallyFrozenAtoms);
1591                     break;
1592                 case (IntegrationAlgorithm::SD1):
1593                     do_update_sd(start_th,
1594                                  end_th,
1595                                  dt,
1596                                  step,
1597                                  x_rvec,
1598                                  xp_rvec,
1599                                  v_rvec,
1600                                  f_rvec,
1601                                  gmx::arrayRefFromArray(inputRecord.opts.nFreeze, inputRecord.opts.ngfrz),
1602                                  invMass,
1603                                  ptype,
1604                                  cFREEZE_,
1605                                  cTC_,
1606                                  inputRecord.ld_seed,
1607                                  cr,
1608                                  sd_,
1609                                  haveConstraints);
1610                     break;
1611                 case (IntegrationAlgorithm::BD):
1612                     do_update_bd(start_th,
1613                                  end_th,
1614                                  dt,
1615                                  step,
1616                                  x_rvec,
1617                                  xp_rvec,
1618                                  v_rvec,
1619                                  f_rvec,
1620                                  gmx::arrayRefFromArray(inputRecord.opts.nFreeze, inputRecord.opts.ngfrz),
1621                                  invMass,
1622                                  ptype,
1623                                  cFREEZE_,
1624                                  cTC_,
1625                                  inputRecord.bd_fric,
1626                                  sd_.bd_rf.data(),
1627                                  inputRecord.ld_seed,
1628                                  DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1629                     break;
1630                 case (IntegrationAlgorithm::VV):
1631                 case (IntegrationAlgorithm::VVAK):
1632                 {
1633                     gmx_bool bExtended = (inputRecord.etc == TemperatureCoupling::NoseHoover
1634                                           || inputRecord.epc == PressureCoupling::ParrinelloRahman
1635                                           || inputRecord.epc == PressureCoupling::Mttk);
1636
1637                     /* assuming barostat coupled to group 0 */
1638                     real alpha = 1.0 + DIM / static_cast<real>(inputRecord.opts.nrdf[0]);
1639                     switch (updatePart)
1640                     {
1641                         case etrtVELOCITY1:
1642                         case etrtVELOCITY2:
1643                             do_update_vv_vel(start_th,
1644                                              end_th,
1645                                              dt,
1646                                              gmx::arrayRefFromArray(inputRecord.opts.nFreeze,
1647                                                                     inputRecord.opts.ngfrz),
1648                                              invMass,
1649                                              ptype,
1650                                              cFREEZE_,
1651                                              v_rvec,
1652                                              f_rvec,
1653                                              bExtended,
1654                                              state->veta,
1655                                              alpha);
1656                             break;
1657                         case etrtPOSITION:
1658                             do_update_vv_pos(start_th,
1659                                              end_th,
1660                                              dt,
1661                                              gmx::arrayRefFromArray(inputRecord.opts.nFreeze,
1662                                                                     inputRecord.opts.ngfrz),
1663                                              ptype,
1664                                              cFREEZE_,
1665                                              x_rvec,
1666                                              xp_rvec,
1667                                              v_rvec,
1668                                              bExtended,
1669                                              state->veta);
1670                             break;
1671                     }
1672                     break;
1673                 }
1674                 default: gmx_fatal(FARGS, "Don't know how to update coordinates");
1675             }
1676         }
1677         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1678     }
1679 }
1680
1681 void Update::Impl::update_for_constraint_virial(const t_inputrec&         inputRecord,
1682                                                 int                       homenr,
1683                                                 bool                      havePartiallyFrozenAtoms,
1684                                                 gmx::ArrayRef<const real> invmass,
1685                                                 gmx::ArrayRef<const rvec> invMassPerDim,
1686                                                 const t_state&            state,
1687                                                 const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
1688                                                 const gmx_ekindata_t& ekind)
1689 {
1690     GMX_ASSERT(inputRecord.eI == IntegrationAlgorithm::MD || inputRecord.eI == IntegrationAlgorithm::SD1,
1691                "Only leap-frog is supported here");
1692
1693     // Cast to real for faster code, no loss in precision
1694     const real dt = inputRecord.delta_t;
1695
1696     const int nth = gmx_omp_nthreads_get(ModuleMultiThread::Update);
1697
1698 #pragma omp parallel for num_threads(nth) schedule(static)
1699     for (int th = 0; th < nth; th++)
1700     {
1701         try
1702         {
1703             int start_th, end_th;
1704             getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
1705
1706             const rvec* x_rvec  = state.x.rvec_array();
1707             rvec*       xp_rvec = xp_.rvec_array();
1708             const rvec* v_rvec  = state.v.rvec_array();
1709             const rvec* f_rvec  = as_rvec_array(f.unpaddedConstArrayRef().data());
1710
1711             doUpdateMDDoNotUpdateVelocities(
1712                     start_th, end_th, dt, x_rvec, xp_rvec, v_rvec, f_rvec, havePartiallyFrozenAtoms, invmass, invMassPerDim, ekind);
1713         }
1714         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1715     }
1716 }