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