src/gromacs/mdlib/update.cpp

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