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