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