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