2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
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"
87 using namespace gmx; // TODO: Remove when this file is moved into gmx namespace
102 std::vector<real> bd_rf;
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;
110 explicit gmx_stochd_t(const t_inputrec& inputRecord);
113 //! pImpled implementation for Update
118 Impl(const t_inputrec& inputRecord, BoxDeformation* boxDeformation);
122 void update_coords(const t_inputrec& inputRecord,
126 const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
128 const gmx_ekindata_t* ekind,
132 bool haveConstraints);
134 void finish_update(const t_inputrec& inputRecord,
137 gmx_wallcycle_t wcycle,
138 bool haveConstraints);
140 void update_sd_second_half(const t_inputrec& inputRecord,
147 gmx_wallcycle_t wcycle,
148 gmx::Constraints* constr,
152 void update_temperature_constants(const t_inputrec& inputRecord);
154 const std::vector<bool>& getAndersenRandomizeGroup() const { return sd_.randomize_group; }
156 const std::vector<real>& getBoltzmanFactor() const { return sd_.boltzfac; }
158 PaddedVector<RVec>* xp() { return &xp_; }
160 BoxDeformation* deform() const { return deform_; }
163 //! stochastic dynamics struct
165 //! xprime for constraint algorithms
166 PaddedVector<RVec> xp_;
167 //! Box deformation handler (or nullptr if inactive).
168 BoxDeformation* deform_ = nullptr;
171 Update::Update(const t_inputrec& inputRecord, BoxDeformation* boxDeformation) :
172 impl_(new Impl(inputRecord, boxDeformation)){};
176 const std::vector<bool>& Update::getAndersenRandomizeGroup() const
178 return impl_->getAndersenRandomizeGroup();
181 const std::vector<real>& Update::getBoltzmanFactor() const
183 return impl_->getBoltzmanFactor();
186 PaddedVector<RVec>* Update::xp()
191 BoxDeformation* Update::deform() const
193 return impl_->deform();
196 void Update::update_coords(const t_inputrec& inputRecord,
200 const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
202 const gmx_ekindata_t* ekind,
206 const bool haveConstraints)
208 return impl_->update_coords(inputRecord, step, md, state, f, fcd, ekind, M, updatePart, cr,
212 void Update::finish_update(const t_inputrec& inputRecord,
215 gmx_wallcycle_t wcycle,
216 const bool haveConstraints)
218 return impl_->finish_update(inputRecord, md, state, wcycle, haveConstraints);
221 void Update::update_sd_second_half(const t_inputrec& inputRecord,
228 gmx_wallcycle_t wcycle,
229 gmx::Constraints* constr,
233 return impl_->update_sd_second_half(inputRecord, step, dvdlambda, md, state, cr, nrnb, wcycle,
234 constr, do_log, do_ene);
237 void Update::update_temperature_constants(const t_inputrec& inputRecord)
239 return impl_->update_temperature_constants(inputRecord);
242 /*! \brief Sets the velocities of virtual sites to zero */
243 static void clearVsiteVelocities(int start, int nrend, const unsigned short* particleType, rvec* gmx_restrict v)
245 for (int a = start; a < nrend; a++)
247 if (particleType[a] == eptVSite)
254 /*! \brief Sets the number of different temperature coupling values */
255 enum class NumTempScaleValues
257 single, //!< Single T-scaling value (either one group or all values =1)
258 multiple //!< Multiple T-scaling values, need to use T-group indices
261 /*! \brief Sets if to apply no or diagonal Parrinello-Rahman pressure scaling
263 * Note that this enum is only used in updateMDLeapfrogSimple(), which does
264 * not handle fully anistropic Parrinello-Rahman scaling, so we only have
265 * options \p no and \p diagonal here and no anistropic option.
267 enum class ApplyParrinelloRahmanVScaling
269 no, //!< Do not apply velocity scaling (not a PR-coupling run or step)
270 diagonal //!< Apply velocity scaling using a diagonal matrix
273 /*! \brief Integrate using leap-frog with T-scaling and optionally diagonal Parrinello-Rahman p-coupling
275 * \tparam numTempScaleValues The number of different T-couple values
276 * \tparam applyPRVScaling Apply Parrinello-Rahman velocity scaling
277 * \param[in] start Index of first atom to update
278 * \param[in] nrend Last atom to update: \p nrend - 1
279 * \param[in] dt The time step
280 * \param[in] dtPressureCouple Time step for pressure coupling
281 * \param[in] invMassPerDim 1/mass per atom and dimension
282 * \param[in] tcstat Temperature coupling information
283 * \param[in] cTC T-coupling group index per atom
284 * \param[in] pRVScaleMatrixDiagonal Parrinello-Rahman v-scale matrix diagonal
285 * \param[in] x Input coordinates
286 * \param[out] xprime Updated coordinates
287 * \param[inout] v Velocities
288 * \param[in] f Forces
290 * We expect this template to get good SIMD acceleration by most compilers,
291 * unlike the more complex general template.
292 * Note that we might get even better SIMD acceleration when we introduce
293 * aligned (and padded) memory, possibly with some hints for the compilers.
295 template<NumTempScaleValues numTempScaleValues, ApplyParrinelloRahmanVScaling applyPRVScaling>
296 static void updateMDLeapfrogSimple(int start,
299 real dtPressureCouple,
300 const rvec* gmx_restrict invMassPerDim,
301 gmx::ArrayRef<const t_grp_tcstat> tcstat,
302 const unsigned short* cTC,
303 const rvec pRVScaleMatrixDiagonal,
304 const rvec* gmx_restrict x,
305 rvec* gmx_restrict xprime,
306 rvec* gmx_restrict v,
307 const rvec* gmx_restrict f)
311 if (numTempScaleValues == NumTempScaleValues::single)
313 lambdaGroup = tcstat[0].lambda;
316 for (int a = start; a < nrend; a++)
318 if (numTempScaleValues == NumTempScaleValues::multiple)
320 lambdaGroup = tcstat[cTC[a]].lambda;
323 for (int d = 0; d < DIM; d++)
325 /* Note that using rvec invMassPerDim results in more efficient
326 * SIMD code, but this increases the cache pressure.
327 * For large systems with PME on the CPU this slows down the
328 * (then already slow) update by 20%. If all data remains in cache,
329 * using rvec is much faster.
331 real vNew = lambdaGroup * v[a][d] + f[a][d] * invMassPerDim[a][d] * dt;
333 if (applyPRVScaling == ApplyParrinelloRahmanVScaling::diagonal)
335 vNew -= dtPressureCouple * pRVScaleMatrixDiagonal[d] * v[a][d];
338 xprime[a][d] = x[a][d] + vNew * dt;
343 #if GMX_SIMD && GMX_SIMD_HAVE_REAL
344 # define GMX_HAVE_SIMD_UPDATE 1
346 # define GMX_HAVE_SIMD_UPDATE 0
349 #if GMX_HAVE_SIMD_UPDATE
351 /*! \brief Load (aligned) the contents of GMX_SIMD_REAL_WIDTH rvec elements sequentially into 3 SIMD registers
353 * The loaded output is:
354 * \p r0: { r[index][XX], r[index][YY], ... }
356 * \p r2: { ..., r[index+GMX_SIMD_REAL_WIDTH-1][YY], r[index+GMX_SIMD_REAL_WIDTH-1][ZZ] }
358 * \param[in] r Real to an rvec array, has to be aligned to SIMD register width
359 * \param[in] index Index of the first rvec triplet of reals to load
360 * \param[out] r0 Pointer to first SIMD register
361 * \param[out] r1 Pointer to second SIMD register
362 * \param[out] r2 Pointer to third SIMD register
364 static inline void simdLoadRvecs(const rvec* r, int index, SimdReal* r0, SimdReal* r1, SimdReal* r2)
366 const real* realPtr = r[index];
368 GMX_ASSERT(isSimdAligned(realPtr), "Pointer should be SIMD aligned");
370 *r0 = simdLoad(realPtr + 0 * GMX_SIMD_REAL_WIDTH);
371 *r1 = simdLoad(realPtr + 1 * GMX_SIMD_REAL_WIDTH);
372 *r2 = simdLoad(realPtr + 2 * GMX_SIMD_REAL_WIDTH);
375 /*! \brief Store (aligned) 3 SIMD registers sequentially to GMX_SIMD_REAL_WIDTH rvec elements
377 * The stored output is:
378 * \p r[index] = { { r0[0], r0[1], ... }
380 * \p r[index+GMX_SIMD_REAL_WIDTH-1] = { ... , r2[GMX_SIMD_REAL_WIDTH-2], r2[GMX_SIMD_REAL_WIDTH-1] }
382 * \param[out] r Pointer to an rvec array, has to be aligned to SIMD register width
383 * \param[in] index Index of the first rvec triplet of reals to store to
384 * \param[in] r0 First SIMD register
385 * \param[in] r1 Second SIMD register
386 * \param[in] r2 Third SIMD register
388 static inline void simdStoreRvecs(rvec* r, int index, SimdReal r0, SimdReal r1, SimdReal r2)
390 real* realPtr = r[index];
392 GMX_ASSERT(isSimdAligned(realPtr), "Pointer should be SIMD aligned");
394 store(realPtr + 0 * GMX_SIMD_REAL_WIDTH, r0);
395 store(realPtr + 1 * GMX_SIMD_REAL_WIDTH, r1);
396 store(realPtr + 2 * GMX_SIMD_REAL_WIDTH, r2);
399 /*! \brief Integrate using leap-frog with single group T-scaling and SIMD
401 * \param[in] start Index of first atom to update
402 * \param[in] nrend Last atom to update: \p nrend - 1
403 * \param[in] dt The time step
404 * \param[in] invMass 1/mass per atom
405 * \param[in] tcstat Temperature coupling information
406 * \param[in] x Input coordinates
407 * \param[out] xprime Updated coordinates
408 * \param[inout] v Velocities
409 * \param[in] f Forces
411 static void updateMDLeapfrogSimpleSimd(int start,
414 const real* gmx_restrict invMass,
415 gmx::ArrayRef<const t_grp_tcstat> tcstat,
416 const rvec* gmx_restrict x,
417 rvec* gmx_restrict xprime,
418 rvec* gmx_restrict v,
419 const rvec* gmx_restrict f)
421 SimdReal timestep(dt);
422 SimdReal lambdaSystem(tcstat[0].lambda);
424 /* We declare variables here, since code is often slower when declaring them inside the loop */
426 /* Note: We should implement a proper PaddedVector, so we don't need this check */
427 GMX_ASSERT(isSimdAligned(invMass), "invMass should be aligned");
429 for (int a = start; a < nrend; a += GMX_SIMD_REAL_WIDTH)
431 SimdReal invMass0, invMass1, invMass2;
432 expandScalarsToTriplets(simdLoad(invMass + a), &invMass0, &invMass1, &invMass2);
436 simdLoadRvecs(v, a, &v0, &v1, &v2);
437 simdLoadRvecs(f, a, &f0, &f1, &f2);
439 v0 = fma(f0 * invMass0, timestep, lambdaSystem * v0);
440 v1 = fma(f1 * invMass1, timestep, lambdaSystem * v1);
441 v2 = fma(f2 * invMass2, timestep, lambdaSystem * v2);
443 simdStoreRvecs(v, a, v0, v1, v2);
446 simdLoadRvecs(x, a, &x0, &x1, &x2);
448 SimdReal xprime0 = fma(v0, timestep, x0);
449 SimdReal xprime1 = fma(v1, timestep, x1);
450 SimdReal xprime2 = fma(v2, timestep, x2);
452 simdStoreRvecs(xprime, a, xprime0, xprime1, xprime2);
456 #endif // GMX_HAVE_SIMD_UPDATE
458 /*! \brief Sets the NEMD acceleration type */
459 enum class AccelerationType
466 /*! \brief Integrate using leap-frog with support for everything.
468 * \tparam accelerationType Type of NEMD acceleration.
469 * \param[in] start Index of first atom to update.
470 * \param[in] nrend Last atom to update: \p nrend - 1.
471 * \param[in] doNoseHoover If to apply Nose-Hoover T-coupling.
472 * \param[in] dt The time step.
473 * \param[in] dtPressureCouple Time step for pressure coupling, is 0 when no pressure
474 * coupling should be applied at this step.
475 * \param[in] accel Acceleration per group.
476 * \param[in] md Atom properties.
477 * \param[in] ekind Kinetic energy data.
478 * \param[in] box The box dimensions.
479 * \param[in] x Input coordinates.
480 * \param[out] xprime Updated coordinates.
481 * \param[inout] v Velocities.
482 * \param[in] f Forces.
483 * \param[in] nh_vxi Nose-Hoover velocity scaling factors.
484 * \param[in] nsttcouple Frequency of the temperature coupling steps.
485 * \param[in] M Parrinello-Rahman scaling matrix.
487 template<AccelerationType accelerationType>
488 static void updateMDLeapfrogGeneral(int start,
492 real dtPressureCouple,
495 const gmx_ekindata_t* ekind,
497 const rvec* gmx_restrict x,
498 rvec* gmx_restrict xprime,
499 rvec* gmx_restrict v,
500 const rvec* gmx_restrict f,
501 const double* gmx_restrict nh_vxi,
502 const int nsttcouple,
505 /* This is a version of the leap-frog integrator that supports
506 * all combinations of T-coupling, P-coupling and NEMD.
507 * Nose-Hoover uses the reversible leap-frog integrator from
508 * Holian et al. Phys Rev E 52(3) : 2338, 1995
511 gmx::ArrayRef<const t_grp_tcstat> tcstat = ekind->tcstat;
512 gmx::ArrayRef<const t_grp_acc> grpstat = ekind->grpstat;
513 const unsigned short* cTC = md->cTC;
514 const unsigned short* cACC = md->cACC;
516 const rvec* gmx_restrict invMassPerDim = md->invMassPerDim;
518 /* Initialize group values, changed later when multiple groups are used */
522 real omega_Z = 2 * static_cast<real>(M_PI) / box[ZZ][ZZ];
524 for (int n = start; n < nrend; n++)
530 real lg = tcstat[gt].lambda;
533 real cosineZ, vCosine;
534 #ifdef __INTEL_COMPILER
535 # pragma warning(disable : 280)
537 switch (accelerationType)
539 case AccelerationType::none: copy_rvec(v[n], vRel); break;
540 case AccelerationType::group:
545 /* Avoid scaling the group velocity */
546 rvec_sub(v[n], grpstat[ga].u, vRel);
548 case AccelerationType::cosine:
549 cosineZ = std::cos(x[n][ZZ] * omega_Z);
550 vCosine = cosineZ * ekind->cosacc.vcos;
551 /* Avoid scaling the cosine profile velocity */
552 copy_rvec(v[n], vRel);
560 /* Here we account for multiple time stepping, by increasing
561 * the Nose-Hoover correction by nsttcouple
562 * TODO: This can be pre-computed.
564 factorNH = 0.5 * nsttcouple * dt * nh_vxi[gt];
567 for (int d = 0; d < DIM; d++)
569 real vNew = (lg * vRel[d]
570 + (f[n][d] * invMassPerDim[n][d] * dt - factorNH * vRel[d]
571 - dtPressureCouple * iprod(M[d], vRel)))
573 switch (accelerationType)
575 case AccelerationType::none: break;
576 case AccelerationType::group:
577 /* Add back the mean velocity and apply acceleration */
578 vNew += grpstat[ga].u[d] + accel[ga][d] * dt;
580 case AccelerationType::cosine:
583 /* Add back the mean velocity and apply acceleration */
584 vNew += vCosine + cosineZ * ekind->cosacc.cos_accel * dt;
589 xprime[n][d] = x[n][d] + vNew * dt;
594 /*! \brief Handles the Leap-frog MD x and v integration */
595 static void do_update_md(int start,
599 const rvec* gmx_restrict x,
600 rvec* gmx_restrict xprime,
601 rvec* gmx_restrict v,
602 const rvec* gmx_restrict f,
603 const rvec* gmx_restrict accel,
606 const int nsttcouple,
607 const int nstpcouple,
609 const gmx_ekindata_t* ekind,
611 const double* gmx_restrict nh_vxi,
614 GMX_ASSERT(nrend == start || xprime != x,
615 "For SIMD optimization certain compilers need to have xprime != x");
617 /* Note: Berendsen pressure scaling is handled after do_update_md() */
618 bool doTempCouple = (etc != etcNO && do_per_step(step + nsttcouple - 1, nsttcouple));
619 bool doNoseHoover = (etc == etcNOSEHOOVER && doTempCouple);
620 bool doParrinelloRahman =
621 (epc == epcPARRINELLORAHMAN && do_per_step(step + nstpcouple - 1, nstpcouple));
622 bool doPROffDiagonal = (doParrinelloRahman && (M[YY][XX] != 0 || M[ZZ][XX] != 0 || M[ZZ][YY] != 0));
624 real dtPressureCouple = (doParrinelloRahman ? nstpcouple * dt : 0);
626 /* NEMD (also cosine) acceleration is applied in updateMDLeapFrogGeneral */
627 bool doAcceleration = (ekind->bNEMD || ekind->cosacc.cos_accel != 0);
629 if (doNoseHoover || doPROffDiagonal || doAcceleration)
632 if (!doParrinelloRahman)
634 /* We should not apply PR scaling at this step */
644 updateMDLeapfrogGeneral<AccelerationType::none>(start, nrend, doNoseHoover, dt,
645 dtPressureCouple, accel, md, ekind, box, x,
646 xprime, v, f, nh_vxi, nsttcouple, stepM);
648 else if (ekind->bNEMD)
650 updateMDLeapfrogGeneral<AccelerationType::group>(
651 start, nrend, doNoseHoover, dt, dtPressureCouple, accel, md, ekind, box, x,
652 xprime, v, f, nh_vxi, nsttcouple, stepM);
656 updateMDLeapfrogGeneral<AccelerationType::cosine>(
657 start, nrend, doNoseHoover, dt, dtPressureCouple, accel, md, ekind, box, x,
658 xprime, v, f, nh_vxi, nsttcouple, stepM);
663 /* Use a simple and thus more efficient integration loop. */
664 /* The simple loop does not check for particle type (so it can
665 * be vectorized), which means we need to clear the velocities
666 * of virtual sites in advance, when present. Note that vsite
667 * velocities are computed after update and constraints from
668 * their displacement.
672 /* Note: The overhead of this loop is completely neligible */
673 clearVsiteVelocities(start, nrend, md->ptype, v);
676 /* We determine if we have a single T-coupling lambda value for all
677 * atoms. That allows for better SIMD acceleration in the template.
678 * If we do not do temperature coupling (in the run or this step),
679 * all scaling values are 1, so we effectively have a single value.
681 bool haveSingleTempScaleValue = (!doTempCouple || ekind->ngtc == 1);
683 /* Extract some pointers needed by all cases */
684 const unsigned short* cTC = md->cTC;
685 gmx::ArrayRef<const t_grp_tcstat> tcstat = ekind->tcstat;
686 const rvec* invMassPerDim = md->invMassPerDim;
688 if (doParrinelloRahman)
690 GMX_ASSERT(!doPROffDiagonal,
691 "updateMDLeapfrogSimple only support diagonal Parrinello-Rahman scaling "
695 for (int d = 0; d < DIM; d++)
700 if (haveSingleTempScaleValue)
702 updateMDLeapfrogSimple<NumTempScaleValues::single, ApplyParrinelloRahmanVScaling::diagonal>(
703 start, nrend, dt, dtPressureCouple, invMassPerDim, tcstat, cTC, diagM, x,
708 updateMDLeapfrogSimple<NumTempScaleValues::multiple, ApplyParrinelloRahmanVScaling::diagonal>(
709 start, nrend, dt, dtPressureCouple, invMassPerDim, tcstat, cTC, diagM, x,
715 if (haveSingleTempScaleValue)
717 /* Note that modern compilers are pretty good at vectorizing
718 * updateMDLeapfrogSimple(). But the SIMD version will still
719 * be faster because invMass lowers the cache pressure
720 * compared to invMassPerDim.
722 #if GMX_HAVE_SIMD_UPDATE
723 /* Check if we can use invmass instead of invMassPerDim */
724 if (!md->havePartiallyFrozenAtoms)
726 updateMDLeapfrogSimpleSimd(start, nrend, dt, md->invmass, tcstat, x, xprime, v, f);
731 updateMDLeapfrogSimple<NumTempScaleValues::single, ApplyParrinelloRahmanVScaling::no>(
732 start, nrend, dt, dtPressureCouple, invMassPerDim, tcstat, cTC, nullptr,
738 updateMDLeapfrogSimple<NumTempScaleValues::multiple, ApplyParrinelloRahmanVScaling::no>(
739 start, nrend, dt, dtPressureCouple, invMassPerDim, tcstat, cTC, nullptr, x,
746 static void do_update_vv_vel(int start,
750 const ivec nFreeze[],
751 const real invmass[],
752 const unsigned short ptype[],
753 const unsigned short cFREEZE[],
754 const unsigned short cACC[],
767 g = 0.25 * dt * veta * alpha;
769 mv2 = gmx::series_sinhx(g);
776 for (n = start; n < nrend; n++)
778 real w_dt = invmass[n] * dt;
788 for (d = 0; d < DIM; d++)
790 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
792 v[n][d] = mv1 * (mv1 * v[n][d] + 0.5 * (w_dt * mv2 * f[n][d])) + 0.5 * accel[ga][d] * dt;
800 } /* do_update_vv_vel */
802 static void do_update_vv_pos(int start,
805 const ivec nFreeze[],
806 const unsigned short ptype[],
807 const unsigned short cFREEZE[],
818 /* Would it make more sense if Parrinello-Rahman was put here? */
823 mr2 = gmx::series_sinhx(g);
831 for (n = start; n < nrend; n++)
839 for (d = 0; d < DIM; d++)
841 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
843 xprime[n][d] = mr1 * (mr1 * x[n][d] + mr2 * dt * v[n][d]);
847 xprime[n][d] = x[n][d];
851 } /* do_update_vv_pos */
853 gmx_stochd_t::gmx_stochd_t(const t_inputrec& inputRecord)
855 const t_grpopts* opts = &inputRecord.opts;
856 const int ngtc = opts->ngtc;
858 if (inputRecord.eI == eiBD)
862 else if (EI_SD(inputRecord.eI))
867 for (int gt = 0; gt < ngtc; gt++)
869 if (opts->tau_t[gt] > 0)
871 sdc[gt].em = std::exp(-inputRecord.delta_t / opts->tau_t[gt]);
875 /* No friction and noise on this group */
880 else if (ETC_ANDERSEN(inputRecord.etc))
882 randomize_group.resize(ngtc);
883 boltzfac.resize(ngtc);
885 /* for now, assume that all groups, if randomized, are randomized at the same rate, i.e. tau_t is the same. */
886 /* since constraint groups don't necessarily match up with temperature groups! This is checked in readir.c */
888 for (int gt = 0; gt < ngtc; gt++)
890 real reft = std::max<real>(0, opts->ref_t[gt]);
891 if ((opts->tau_t[gt] > 0)
892 && (reft > 0)) /* tau_t or ref_t = 0 means that no randomization is done */
894 randomize_group[gt] = true;
895 boltzfac[gt] = BOLTZ * opts->ref_t[gt];
899 randomize_group[gt] = false;
905 void Update::Impl::update_temperature_constants(const t_inputrec& inputRecord)
907 if (inputRecord.eI == eiBD)
909 if (inputRecord.bd_fric != 0)
911 for (int gt = 0; gt < inputRecord.opts.ngtc; gt++)
913 sd_.bd_rf[gt] = std::sqrt(2.0 * BOLTZ * inputRecord.opts.ref_t[gt]
914 / (inputRecord.bd_fric * inputRecord.delta_t));
919 for (int gt = 0; gt < inputRecord.opts.ngtc; gt++)
921 sd_.bd_rf[gt] = std::sqrt(2.0 * BOLTZ * inputRecord.opts.ref_t[gt]);
925 if (inputRecord.eI == eiSD1)
927 for (int gt = 0; gt < inputRecord.opts.ngtc; gt++)
929 real kT = BOLTZ * inputRecord.opts.ref_t[gt];
930 /* The mass is accounted for later, since this differs per atom */
931 sd_.sdsig[gt].V = std::sqrt(kT * (1 - sd_.sdc[gt].em * sd_.sdc[gt].em));
936 Update::Impl::Impl(const t_inputrec& inputRecord, BoxDeformation* boxDeformation) :
938 deform_(boxDeformation)
940 update_temperature_constants(inputRecord);
941 xp_.resizeWithPadding(0);
944 void Update::setNumAtoms(int numAtoms)
947 impl_->xp()->resizeWithPadding(numAtoms);
950 /*! \brief Sets the SD update type */
951 enum class SDUpdate : int
954 FrictionAndNoiseOnly,
958 /*! \brief SD integrator update
960 * Two phases are required in the general case of a constrained
961 * update, the first phase from the contribution of forces, before
962 * applying constraints, and then a second phase applying the friction
963 * and noise, and then further constraining. For details, see
966 * Without constraints, the two phases can be combined, for
969 * Thus three instantiations of this templated function will be made,
970 * two with only one contribution, and one with both contributions. */
971 template<SDUpdate updateType>
972 static void doSDUpdateGeneral(const gmx_stochd_t& sd,
977 const ivec nFreeze[],
978 const real invmass[],
979 const unsigned short ptype[],
980 const unsigned short cFREEZE[],
981 const unsigned short cACC[],
982 const unsigned short cTC[],
991 // cTC, cACC and cFREEZE can be nullptr any time, but various
992 // instantiations do not make sense with particular pointer
994 if (updateType == SDUpdate::ForcesOnly)
996 GMX_ASSERT(f != nullptr, "SD update with only forces requires forces");
997 GMX_ASSERT(cTC == nullptr, "SD update with only forces cannot handle temperature groups");
999 if (updateType == SDUpdate::FrictionAndNoiseOnly)
1001 GMX_ASSERT(f == nullptr, "SD update with only noise cannot handle forces");
1002 GMX_ASSERT(cACC == nullptr, "SD update with only noise cannot handle acceleration groups");
1004 if (updateType == SDUpdate::Combined)
1006 GMX_ASSERT(f != nullptr, "SD update with forces and noise requires forces");
1009 // Even 0 bits internal counter gives 2x64 ints (more than enough for three table lookups)
1010 gmx::ThreeFry2x64<0> rng(seed, gmx::RandomDomain::UpdateCoordinates);
1011 gmx::TabulatedNormalDistribution<real, 14> dist;
1013 for (int n = start; n < nrend; n++)
1015 int globalAtomIndex = gatindex ? gatindex[n] : n;
1016 rng.restart(step, globalAtomIndex);
1019 real inverseMass = invmass[n];
1020 real invsqrtMass = std::sqrt(inverseMass);
1022 int freezeGroup = cFREEZE ? cFREEZE[n] : 0;
1023 int accelerationGroup = cACC ? cACC[n] : 0;
1024 int temperatureGroup = cTC ? cTC[n] : 0;
1026 for (int d = 0; d < DIM; d++)
1028 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[freezeGroup][d])
1030 if (updateType == SDUpdate::ForcesOnly)
1032 real vn = v[n][d] + (inverseMass * f[n][d] + accel[accelerationGroup][d]) * dt;
1034 // Simple position update.
1035 xprime[n][d] = x[n][d] + v[n][d] * dt;
1037 else if (updateType == SDUpdate::FrictionAndNoiseOnly)
1040 v[n][d] = (vn * sd.sdc[temperatureGroup].em
1041 + invsqrtMass * sd.sdsig[temperatureGroup].V * dist(rng));
1042 // The previous phase already updated the
1043 // positions with a full v*dt term that must
1044 // now be half removed.
1045 xprime[n][d] = xprime[n][d] + 0.5 * (v[n][d] - vn) * dt;
1049 real vn = v[n][d] + (inverseMass * f[n][d] + accel[accelerationGroup][d]) * dt;
1050 v[n][d] = (vn * sd.sdc[temperatureGroup].em
1051 + invsqrtMass * sd.sdsig[temperatureGroup].V * dist(rng));
1052 // Here we include half of the friction+noise
1053 // update of v into the position update.
1054 xprime[n][d] = x[n][d] + 0.5 * (vn + v[n][d]) * dt;
1059 // When using constraints, the update is split into
1060 // two phases, but we only need to zero the update of
1061 // virtual, shell or frozen particles in at most one
1063 if (updateType != SDUpdate::FrictionAndNoiseOnly)
1066 xprime[n][d] = x[n][d];
1073 static void do_update_sd(int start,
1077 const rvec* gmx_restrict x,
1078 rvec* gmx_restrict xprime,
1079 rvec* gmx_restrict v,
1080 const rvec* gmx_restrict f,
1082 const ivec nFreeze[],
1083 const real invmass[],
1084 const unsigned short ptype[],
1085 const unsigned short cFREEZE[],
1086 const unsigned short cACC[],
1087 const unsigned short cTC[],
1089 const t_commrec* cr,
1090 const gmx_stochd_t& sd,
1091 bool haveConstraints)
1093 if (haveConstraints)
1095 // With constraints, the SD update is done in 2 parts
1096 doSDUpdateGeneral<SDUpdate::ForcesOnly>(sd, start, nrend, dt, accel, nFreeze, invmass,
1097 ptype, cFREEZE, cACC, nullptr, x, xprime, v, f,
1098 step, seed, nullptr);
1102 doSDUpdateGeneral<SDUpdate::Combined>(
1103 sd, start, nrend, dt, accel, nFreeze, invmass, ptype, cFREEZE, cACC, cTC, x, xprime,
1104 v, f, step, seed, DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1108 static void do_update_bd(int start,
1112 const rvec* gmx_restrict x,
1113 rvec* gmx_restrict xprime,
1114 rvec* gmx_restrict v,
1115 const rvec* gmx_restrict f,
1116 const ivec nFreeze[],
1117 const real invmass[],
1118 const unsigned short ptype[],
1119 const unsigned short cFREEZE[],
1120 const unsigned short cTC[],
1121 real friction_coefficient,
1124 const int* gatindex)
1126 /* note -- these appear to be full step velocities . . . */
1131 // Use 1 bit of internal counters to give us 2*2 64-bits values per stream
1132 // Each 64-bit value is enough for 4 normal distribution table numbers.
1133 gmx::ThreeFry2x64<0> rng(seed, gmx::RandomDomain::UpdateCoordinates);
1134 gmx::TabulatedNormalDistribution<real, 14> dist;
1136 if (friction_coefficient != 0)
1138 invfr = 1.0 / friction_coefficient;
1141 for (n = start; (n < nrend); n++)
1143 int ng = gatindex ? gatindex[n] : n;
1145 rng.restart(step, ng);
1156 for (d = 0; (d < DIM); d++)
1158 if ((ptype[n] != eptVSite) && (ptype[n] != eptShell) && !nFreeze[gf][d])
1160 if (friction_coefficient != 0)
1162 vn = invfr * f[n][d] + rf[gt] * dist(rng);
1166 /* NOTE: invmass = 2/(mass*friction_constant*dt) */
1167 vn = 0.5 * invmass[n] * f[n][d] * dt
1168 + std::sqrt(0.5 * invmass[n]) * rf[gt] * dist(rng);
1172 xprime[n][d] = x[n][d] + vn * dt;
1177 xprime[n][d] = x[n][d];
1183 extern void init_ekinstate(ekinstate_t* ekinstate, const t_inputrec* ir)
1185 ekinstate->ekin_n = ir->opts.ngtc;
1186 snew(ekinstate->ekinh, ekinstate->ekin_n);
1187 snew(ekinstate->ekinf, ekinstate->ekin_n);
1188 snew(ekinstate->ekinh_old, ekinstate->ekin_n);
1189 ekinstate->ekinscalef_nhc.resize(ekinstate->ekin_n);
1190 ekinstate->ekinscaleh_nhc.resize(ekinstate->ekin_n);
1191 ekinstate->vscale_nhc.resize(ekinstate->ekin_n);
1192 ekinstate->dekindl = 0;
1193 ekinstate->mvcos = 0;
1194 ekinstate->hasReadEkinState = false;
1197 void update_ekinstate(ekinstate_t* ekinstate, const gmx_ekindata_t* ekind)
1201 for (i = 0; i < ekinstate->ekin_n; i++)
1203 copy_mat(ekind->tcstat[i].ekinh, ekinstate->ekinh[i]);
1204 copy_mat(ekind->tcstat[i].ekinf, ekinstate->ekinf[i]);
1205 copy_mat(ekind->tcstat[i].ekinh_old, ekinstate->ekinh_old[i]);
1206 ekinstate->ekinscalef_nhc[i] = ekind->tcstat[i].ekinscalef_nhc;
1207 ekinstate->ekinscaleh_nhc[i] = ekind->tcstat[i].ekinscaleh_nhc;
1208 ekinstate->vscale_nhc[i] = ekind->tcstat[i].vscale_nhc;
1211 copy_mat(ekind->ekin, ekinstate->ekin_total);
1212 ekinstate->dekindl = ekind->dekindl;
1213 ekinstate->mvcos = ekind->cosacc.mvcos;
1216 void restore_ekinstate_from_state(const t_commrec* cr, gmx_ekindata_t* ekind, const ekinstate_t* ekinstate)
1222 for (i = 0; i < ekinstate->ekin_n; i++)
1224 copy_mat(ekinstate->ekinh[i], ekind->tcstat[i].ekinh);
1225 copy_mat(ekinstate->ekinf[i], ekind->tcstat[i].ekinf);
1226 copy_mat(ekinstate->ekinh_old[i], ekind->tcstat[i].ekinh_old);
1227 ekind->tcstat[i].ekinscalef_nhc = ekinstate->ekinscalef_nhc[i];
1228 ekind->tcstat[i].ekinscaleh_nhc = ekinstate->ekinscaleh_nhc[i];
1229 ekind->tcstat[i].vscale_nhc = ekinstate->vscale_nhc[i];
1232 copy_mat(ekinstate->ekin_total, ekind->ekin);
1234 ekind->dekindl = ekinstate->dekindl;
1235 ekind->cosacc.mvcos = ekinstate->mvcos;
1236 n = ekinstate->ekin_n;
1241 gmx_bcast(sizeof(n), &n, cr->mpi_comm_mygroup);
1242 for (i = 0; i < n; i++)
1244 gmx_bcast(DIM * DIM * sizeof(ekind->tcstat[i].ekinh[0][0]), ekind->tcstat[i].ekinh[0],
1245 cr->mpi_comm_mygroup);
1246 gmx_bcast(DIM * DIM * sizeof(ekind->tcstat[i].ekinf[0][0]), ekind->tcstat[i].ekinf[0],
1247 cr->mpi_comm_mygroup);
1248 gmx_bcast(DIM * DIM * sizeof(ekind->tcstat[i].ekinh_old[0][0]),
1249 ekind->tcstat[i].ekinh_old[0], cr->mpi_comm_mygroup);
1251 gmx_bcast(sizeof(ekind->tcstat[i].ekinscalef_nhc), &(ekind->tcstat[i].ekinscalef_nhc),
1252 cr->mpi_comm_mygroup);
1253 gmx_bcast(sizeof(ekind->tcstat[i].ekinscaleh_nhc), &(ekind->tcstat[i].ekinscaleh_nhc),
1254 cr->mpi_comm_mygroup);
1255 gmx_bcast(sizeof(ekind->tcstat[i].vscale_nhc), &(ekind->tcstat[i].vscale_nhc),
1256 cr->mpi_comm_mygroup);
1258 gmx_bcast(DIM * DIM * sizeof(ekind->ekin[0][0]), ekind->ekin[0], cr->mpi_comm_mygroup);
1260 gmx_bcast(sizeof(ekind->dekindl), &ekind->dekindl, cr->mpi_comm_mygroup);
1261 gmx_bcast(sizeof(ekind->cosacc.mvcos), &ekind->cosacc.mvcos, cr->mpi_comm_mygroup);
1265 void update_tcouple(int64_t step,
1266 const t_inputrec* inputrec,
1268 gmx_ekindata_t* ekind,
1269 const t_extmass* MassQ,
1270 const t_mdatoms* md)
1273 // This condition was explicitly checked in previous version, but should have never been satisfied
1274 GMX_ASSERT(!(EI_VV(inputrec->eI)
1275 && (inputrecNvtTrotter(inputrec) || inputrecNptTrotter(inputrec)
1276 || inputrecNphTrotter(inputrec))),
1277 "Temperature coupling was requested with velocity verlet and trotter");
1279 bool doTemperatureCoupling = false;
1281 // For VV temperature coupling parameters are updated on the current
1282 // step, for the others - one step before.
1283 if (inputrec->etc == etcNO)
1285 doTemperatureCoupling = false;
1287 else if (EI_VV(inputrec->eI))
1289 doTemperatureCoupling = do_per_step(step, inputrec->nsttcouple);
1293 doTemperatureCoupling = do_per_step(step + inputrec->nsttcouple - 1, inputrec->nsttcouple);
1296 if (doTemperatureCoupling)
1298 real dttc = inputrec->nsttcouple * inputrec->delta_t;
1300 // TODO: berendsen_tcoupl(...), nosehoover_tcoupl(...) and vrescale_tcoupl(...) update
1301 // temperature coupling parameters, which should be reflected in the name of these
1303 switch (inputrec->etc)
1307 berendsen_tcoupl(inputrec, ekind, dttc, state->therm_integral);
1310 nosehoover_tcoupl(&(inputrec->opts), ekind, dttc, state->nosehoover_xi.data(),
1311 state->nosehoover_vxi.data(), MassQ);
1314 vrescale_tcoupl(inputrec, step, ekind, dttc, state->therm_integral.data());
1317 /* rescale in place here */
1318 if (EI_VV(inputrec->eI))
1320 rescale_velocities(ekind, md, 0, md->homenr, state->v.rvec_array());
1325 // Set the T scaling lambda to 1 to have no scaling
1326 // TODO: Do we have to do it on every non-t-couple step?
1327 for (int i = 0; (i < inputrec->opts.ngtc); i++)
1329 ekind->tcstat[i].lambda = 1.0;
1334 void getThreadAtomRange(int numThreads, int threadIndex, int numAtoms, int* startAtom, int* endAtom)
1336 #if GMX_HAVE_SIMD_UPDATE
1337 constexpr int blockSize = GMX_SIMD_REAL_WIDTH;
1339 constexpr int blockSize = 1;
1341 int numBlocks = (numAtoms + blockSize - 1) / blockSize;
1343 *startAtom = ((numBlocks * threadIndex) / numThreads) * blockSize;
1344 *endAtom = ((numBlocks * (threadIndex + 1)) / numThreads) * blockSize;
1345 if (threadIndex == numThreads - 1)
1347 *endAtom = numAtoms;
1351 void update_pcouple_before_coordinates(FILE* fplog,
1353 const t_inputrec* inputrec,
1355 matrix parrinellorahmanMu,
1359 /* Berendsen P-coupling is completely handled after the coordinate update.
1360 * Trotter P-coupling is handled by separate calls to trotter_update().
1362 if (inputrec->epc == epcPARRINELLORAHMAN
1363 && do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
1365 real dtpc = inputrec->nstpcouple * inputrec->delta_t;
1367 parrinellorahman_pcoupl(fplog, step, inputrec, dtpc, state->pres_prev, state->box,
1368 state->box_rel, state->boxv, M, parrinellorahmanMu, bInitStep);
1372 void Update::Impl::update_sd_second_half(const t_inputrec& inputRecord,
1375 const t_mdatoms* md,
1377 const t_commrec* cr,
1379 gmx_wallcycle_t wcycle,
1380 gmx::Constraints* constr,
1388 if (inputRecord.eI == eiSD1)
1390 int homenr = md->homenr;
1392 /* Cast delta_t from double to real to make the integrators faster.
1393 * The only reason for having delta_t double is to get accurate values
1394 * for t=delta_t*step when step is larger than float precision.
1395 * For integration dt the accuracy of real suffices, since with
1396 * integral += dt*integrand the increment is nearly always (much) smaller
1397 * than the integral (and the integrand has real precision).
1399 real dt = inputRecord.delta_t;
1401 wallcycle_start(wcycle, ewcUPDATE);
1403 int nth = gmx_omp_nthreads_get(emntUpdate);
1405 #pragma omp parallel for num_threads(nth) schedule(static)
1406 for (int th = 0; th < nth; th++)
1410 int start_th, end_th;
1411 getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
1413 doSDUpdateGeneral<SDUpdate::FrictionAndNoiseOnly>(
1414 sd_, start_th, end_th, dt, inputRecord.opts.acc, inputRecord.opts.nFreeze,
1415 md->invmass, md->ptype, md->cFREEZE, nullptr, md->cTC, state->x.rvec_array(),
1416 xp_.rvec_array(), state->v.rvec_array(), nullptr, step, inputRecord.ld_seed,
1417 DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1419 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1421 inc_nrnb(nrnb, eNR_UPDATE, homenr);
1422 wallcycle_stop(wcycle, ewcUPDATE);
1424 /* Constrain the coordinates upd->xp for half a time step */
1425 bool computeVirial = false;
1426 constr->apply(do_log, do_ene, step, 1, 0.5, state->x.arrayRefWithPadding(),
1427 xp_.arrayRefWithPadding(), ArrayRef<RVec>(), state->box,
1428 state->lambda[efptBONDED], dvdlambda, state->v.arrayRefWithPadding(),
1429 computeVirial, nullptr, ConstraintVariable::Positions);
1433 void Update::Impl::finish_update(const t_inputrec& inputRecord,
1434 const t_mdatoms* md,
1436 gmx_wallcycle_t wcycle,
1437 const bool haveConstraints)
1439 /* NOTE: Currently we always integrate to a temporary buffer and
1440 * then copy the results back here.
1443 wallcycle_start_nocount(wcycle, ewcUPDATE);
1445 const int homenr = md->homenr;
1446 auto xp = makeConstArrayRef(xp_).subArray(0, homenr);
1447 auto x = makeArrayRef(state->x).subArray(0, homenr);
1449 if (md->havePartiallyFrozenAtoms && haveConstraints)
1451 /* We have atoms that are frozen along some, but not all dimensions,
1452 * then constraints will have moved them also along the frozen dimensions.
1453 * To freeze such degrees of freedom we do not copy them back here.
1455 const ivec* nFreeze = inputRecord.opts.nFreeze;
1457 for (int i = 0; i < homenr; i++)
1459 const int g = md->cFREEZE[i];
1461 for (int d = 0; d < DIM; d++)
1463 if (nFreeze[g][d] == 0)
1472 /* We have no frozen atoms or fully frozen atoms which have not
1473 * been moved by the update, so we can simply copy all coordinates.
1475 int gmx_unused nth = gmx_omp_nthreads_get(emntUpdate);
1476 #pragma omp parallel for num_threads(nth) schedule(static)
1477 for (int i = 0; i < homenr; i++)
1479 // Trivial statement, does not throw
1484 wallcycle_stop(wcycle, ewcUPDATE);
1487 void update_pcouple_after_coordinates(FILE* fplog,
1489 const t_inputrec* inputrec,
1490 const t_mdatoms* md,
1491 const matrix pressure,
1492 const matrix forceVirial,
1493 const matrix constraintVirial,
1494 matrix pressureCouplingMu,
1497 gmx::BoxDeformation* boxDeformation,
1498 const bool scaleCoordinates)
1501 int homenr = md->homenr;
1503 /* Cast to real for faster code, no loss in precision (see comment above) */
1504 real dt = inputrec->delta_t;
1507 /* now update boxes */
1508 switch (inputrec->epc)
1510 case (epcNO): break;
1511 case (epcBERENDSEN):
1512 if (do_per_step(step, inputrec->nstpcouple))
1514 real dtpc = inputrec->nstpcouple * dt;
1515 berendsen_pcoupl(fplog, step, inputrec, dtpc, pressure, state->box, forceVirial,
1516 constraintVirial, pressureCouplingMu, &state->baros_integral);
1517 berendsen_pscale(inputrec, pressureCouplingMu, state->box, state->box_rel, start,
1518 homenr, state->x.rvec_array(), md->cFREEZE, nrnb, scaleCoordinates);
1521 case (epcPARRINELLORAHMAN):
1522 if (do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
1524 /* The box velocities were updated in do_pr_pcoupl,
1525 * but we dont change the box vectors until we get here
1526 * since we need to be able to shift/unshift above.
1528 real dtpc = inputrec->nstpcouple * dt;
1529 for (int i = 0; i < DIM; i++)
1531 for (int m = 0; m <= i; m++)
1533 state->box[i][m] += dtpc * state->boxv[i][m];
1536 preserve_box_shape(inputrec, state->box_rel, state->box);
1538 /* Scale the coordinates */
1539 if (scaleCoordinates)
1541 auto x = state->x.rvec_array();
1542 for (int n = start; n < start + homenr; n++)
1544 tmvmul_ur0(pressureCouplingMu, x[n], x[n]);
1550 switch (inputrec->epct)
1552 case (epctISOTROPIC):
1553 /* DIM * eta = ln V. so DIM*eta_new = DIM*eta_old + DIM*dt*veta =>
1554 ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) =>
1555 Side length scales as exp(veta*dt) */
1557 msmul(state->box, std::exp(state->veta * dt), state->box);
1559 /* Relate veta to boxv. veta = d(eta)/dT = (1/DIM)*1/V dV/dT.
1560 o If we assume isotropic scaling, and box length scaling
1561 factor L, then V = L^DIM (det(M)). So dV/dt = DIM
1562 L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt. The
1563 determinant of B is L^DIM det(M), and the determinant
1564 of dB/dt is (dL/dT)^DIM det (M). veta will be
1565 (det(dB/dT)/det(B))^(1/3). Then since M =
1566 B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */
1568 msmul(state->box, state->veta, state->boxv);
1578 auto localX = makeArrayRef(state->x).subArray(start, homenr);
1579 boxDeformation->apply(localX, state->box, step);
1583 void Update::Impl::update_coords(const t_inputrec& inputRecord,
1585 const t_mdatoms* md,
1587 const gmx::ArrayRefWithPadding<const gmx::RVec>& f,
1588 const t_fcdata* fcd,
1589 const gmx_ekindata_t* ekind,
1592 const t_commrec* cr,
1593 const bool haveConstraints)
1595 /* Running the velocity half does nothing except for velocity verlet */
1596 if ((updatePart == etrtVELOCITY1 || updatePart == etrtVELOCITY2) && !EI_VV(inputRecord.eI))
1598 gmx_incons("update_coords called for velocity without VV integrator");
1601 int homenr = md->homenr;
1603 /* Cast to real for faster code, no loss in precision (see comment above) */
1604 real dt = inputRecord.delta_t;
1606 /* We need to update the NMR restraint history when time averaging is used */
1607 if (state->flags & (1 << estDISRE_RM3TAV))
1609 update_disres_history(fcd, &state->hist);
1611 if (state->flags & (1 << estORIRE_DTAV))
1613 update_orires_history(fcd, &state->hist);
1616 /* ############# START The update of velocities and positions ######### */
1617 int nth = gmx_omp_nthreads_get(emntUpdate);
1619 #pragma omp parallel for num_threads(nth) schedule(static)
1620 for (int th = 0; th < nth; th++)
1624 int start_th, end_th;
1625 getThreadAtomRange(nth, th, homenr, &start_th, &end_th);
1627 const rvec* x_rvec = state->x.rvec_array();
1628 rvec* xp_rvec = xp_.rvec_array();
1629 rvec* v_rvec = state->v.rvec_array();
1630 const rvec* f_rvec = as_rvec_array(f.unpaddedConstArrayRef().data());
1632 switch (inputRecord.eI)
1635 do_update_md(start_th, end_th, dt, step, x_rvec, xp_rvec, v_rvec, f_rvec,
1636 inputRecord.opts.acc, inputRecord.etc, inputRecord.epc,
1637 inputRecord.nsttcouple, inputRecord.nstpcouple, md, ekind,
1638 state->box, state->nosehoover_vxi.data(), M);
1641 do_update_sd(start_th, end_th, dt, step, x_rvec, xp_rvec, v_rvec, f_rvec,
1642 inputRecord.opts.acc, inputRecord.opts.nFreeze, md->invmass,
1643 md->ptype, md->cFREEZE, md->cACC, md->cTC, inputRecord.ld_seed, cr,
1644 sd_, haveConstraints);
1647 do_update_bd(start_th, end_th, dt, step, x_rvec, xp_rvec, v_rvec, f_rvec,
1648 inputRecord.opts.nFreeze, md->invmass, md->ptype, md->cFREEZE,
1649 md->cTC, inputRecord.bd_fric, sd_.bd_rf.data(), inputRecord.ld_seed,
1650 DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1655 gmx_bool bExtended = (inputRecord.etc == etcNOSEHOOVER || inputRecord.epc == epcPARRINELLORAHMAN
1656 || inputRecord.epc == epcMTTK);
1658 /* assuming barostat coupled to group 0 */
1659 real alpha = 1.0 + DIM / static_cast<real>(inputRecord.opts.nrdf[0]);
1664 do_update_vv_vel(start_th, end_th, dt, inputRecord.opts.acc,
1665 inputRecord.opts.nFreeze, md->invmass, md->ptype, md->cFREEZE,
1666 md->cACC, v_rvec, f_rvec, bExtended, state->veta, alpha);
1669 do_update_vv_pos(start_th, end_th, dt, inputRecord.opts.nFreeze,
1670 md->ptype, md->cFREEZE, x_rvec, xp_rvec, v_rvec,
1671 bExtended, state->veta);
1676 default: gmx_fatal(FARGS, "Don't know how to update coordinates");
1679 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
1683 extern gmx_bool update_randomize_velocities(const t_inputrec* ir,
1685 const t_commrec* cr,
1686 const t_mdatoms* md,
1687 gmx::ArrayRef<gmx::RVec> v,
1689 const gmx::Constraints* constr)
1692 real rate = (ir->delta_t) / ir->opts.tau_t[0];
1694 if (ir->etc == etcANDERSEN && constr != nullptr)
1696 /* Currently, Andersen thermostat does not support constrained
1697 systems. Functionality exists in the andersen_tcoupl
1698 function in GROMACS 4.5.7 to allow this combination. That
1699 code could be ported to the current random-number
1700 generation approach, but has not yet been done because of
1701 lack of time and resources. */
1703 "Normal Andersen is currently not supported with constraints, use massive "
1704 "Andersen instead");
1707 /* proceed with andersen if 1) it's fixed probability per
1708 particle andersen or 2) it's massive andersen and it's tau_t/dt */
1709 if ((ir->etc == etcANDERSEN) || do_per_step(step, roundToInt(1.0 / rate)))
1711 andersen_tcoupl(ir, step, cr, md, v, rate, upd->getAndersenRandomizeGroup(),
1712 upd->getBoltzmanFactor());