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,2018, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
47 #include "gromacs/domdec/domdec_struct.h"
48 #include "gromacs/fileio/confio.h"
49 #include "gromacs/gmxlib/network.h"
50 #include "gromacs/gmxlib/nrnb.h"
51 #include "gromacs/listed-forces/disre.h"
52 #include "gromacs/listed-forces/orires.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/invertmatrix.h"
55 #include "gromacs/math/paddedvector.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/math/vecdump.h"
59 #include "gromacs/mdlib/constr.h"
60 #include "gromacs/mdlib/force.h"
61 #include "gromacs/mdlib/gmx_omp_nthreads.h"
62 #include "gromacs/mdlib/mdrun.h"
63 #include "gromacs/mdlib/sim_util.h"
64 #include "gromacs/mdlib/tgroup.h"
65 #include "gromacs/mdtypes/commrec.h"
66 #include "gromacs/mdtypes/group.h"
67 #include "gromacs/mdtypes/inputrec.h"
68 #include "gromacs/mdtypes/md_enums.h"
69 #include "gromacs/mdtypes/state.h"
70 #include "gromacs/pbcutil/boxutilities.h"
71 #include "gromacs/pbcutil/mshift.h"
72 #include "gromacs/pbcutil/pbc.h"
73 #include "gromacs/pulling/pull.h"
74 #include "gromacs/random/tabulatednormaldistribution.h"
75 #include "gromacs/random/threefry.h"
76 #include "gromacs/simd/simd.h"
77 #include "gromacs/timing/wallcycle.h"
78 #include "gromacs/topology/atoms.h"
79 #include "gromacs/utility/exceptions.h"
80 #include "gromacs/utility/fatalerror.h"
81 #include "gromacs/utility/futil.h"
82 #include "gromacs/utility/gmxassert.h"
83 #include "gromacs/utility/gmxomp.h"
84 #include "gromacs/utility/smalloc.h"
86 using namespace gmx
; // TODO: Remove when this file is moved into gmx namespace
101 gmx_sd_sigma_t
*sdsig
;
102 /* andersen temperature control stuff */
103 gmx_bool
*randomize_group
;
110 /* xprime for constraint algorithms */
113 /* Variables for the deform algorithm */
114 gmx_int64_t deformref_step
;
115 matrix deformref_box
;
118 static bool isTemperatureCouplingStep(gmx_int64_t step
, const t_inputrec
*ir
)
120 /* We should only couple after a step where energies were determined (for leapfrog versions)
121 or the step energies are determined, for velocity verlet versions */
131 return ir
->etc
!= etcNO
&&
132 (ir
->nsttcouple
== 1 ||
133 do_per_step(step
+ ir
->nsttcouple
- offset
, ir
->nsttcouple
));
136 static bool isPressureCouplingStep(gmx_int64_t step
, const t_inputrec
*ir
)
138 GMX_ASSERT(ir
->epc
!= epcMTTK
, "MTTK pressure coupling is not handled here");
141 if (ir
->epc
== epcBERENDSEN
)
149 /* We should only couple after a step where pressures were determined */
150 return ir
->epc
!= etcNO
&&
151 (ir
->nstpcouple
== 1 ||
152 do_per_step(step
+ ir
->nstpcouple
- offset
, ir
->nstpcouple
));
155 /*! \brief Sets the velocities of virtual sites to zero */
156 static void clearVsiteVelocities(int start
,
158 const unsigned short *particleType
,
159 rvec
* gmx_restrict v
)
161 for (int a
= start
; a
< nrend
; a
++)
163 if (particleType
[a
] == eptVSite
)
170 /*! \brief Sets the number of different temperature coupling values */
171 enum class NumTempScaleValues
173 single
, //!< Single T-scaling value (either one group or all values =1)
174 multiple
//!< Multiple T-scaling values, need to use T-group indices
177 /*! \brief Sets if to apply no or diagonal Parrinello-Rahman pressure scaling
179 * Note that this enum is only used in updateMDLeapfrogSimple(), which does
180 * not handle fully anistropic Parrinello-Rahman scaling, so we only have
181 * options \p no and \p diagonal here and no anistropic option.
183 enum class ApplyParrinelloRahmanVScaling
185 no
, //!< Do not apply velocity scaling (not a PR-coupling run or step)
186 diagonal
//!< Apply velocity scaling using a diagonal matrix
189 /*! \brief Integrate using leap-frog with T-scaling and optionally diagonal Parrinello-Rahman p-coupling
191 * \tparam numTempScaleValues The number of different T-couple values
192 * \tparam applyPRVScaling Apply Parrinello-Rahman velocity scaling
193 * \param[in] start Index of first atom to update
194 * \param[in] nrend Last atom to update: \p nrend - 1
195 * \param[in] dt The time step
196 * \param[in] dtPressureCouple Time step for pressure coupling
197 * \param[in] invMassPerDim 1/mass per atom and dimension
198 * \param[in] tcstat Temperature coupling information
199 * \param[in] cTC T-coupling group index per atom
200 * \param[in] pRVScaleMatrixDiagonal Parrinello-Rahman v-scale matrix diagonal
201 * \param[in] x Input coordinates
202 * \param[out] xprime Updated coordinates
203 * \param[inout] v Velocities
204 * \param[in] f Forces
206 * We expect this template to get good SIMD acceleration by most compilers,
207 * unlike the more complex general template.
208 * Note that we might get even better SIMD acceleration when we introduce
209 * aligned (and padded) memory, possibly with some hints for the compilers.
211 template<NumTempScaleValues numTempScaleValues
,
212 ApplyParrinelloRahmanVScaling applyPRVScaling
>
214 updateMDLeapfrogSimple(int start
,
217 real dtPressureCouple
,
218 const rvec
* gmx_restrict invMassPerDim
,
219 const t_grp_tcstat
* tcstat
,
220 const unsigned short * cTC
,
221 const rvec pRVScaleMatrixDiagonal
,
222 const rvec
* gmx_restrict x
,
223 rvec
* gmx_restrict xprime
,
224 rvec
* gmx_restrict v
,
225 const rvec
* gmx_restrict f
)
229 if (numTempScaleValues
== NumTempScaleValues::single
)
231 lambdaGroup
= tcstat
[0].lambda
;
234 for (int a
= start
; a
< nrend
; a
++)
236 if (numTempScaleValues
== NumTempScaleValues::multiple
)
238 lambdaGroup
= tcstat
[cTC
[a
]].lambda
;
241 for (int d
= 0; d
< DIM
; d
++)
243 /* Note that using rvec invMassPerDim results in more efficient
244 * SIMD code, but this increases the cache pressure.
245 * For large systems with PME on the CPU this slows down the
246 * (then already slow) update by 20%. If all data remains in cache,
247 * using rvec is much faster.
249 real vNew
= lambdaGroup
*v
[a
][d
] + f
[a
][d
]*invMassPerDim
[a
][d
]*dt
;
251 if (applyPRVScaling
== ApplyParrinelloRahmanVScaling::diagonal
)
253 vNew
-= dtPressureCouple
*pRVScaleMatrixDiagonal
[d
]*v
[a
][d
];
256 xprime
[a
][d
] = x
[a
][d
] + vNew
*dt
;
261 #if GMX_SIMD && GMX_SIMD_HAVE_REAL
262 #define GMX_HAVE_SIMD_UPDATE 1
264 #define GMX_HAVE_SIMD_UPDATE 0
267 #if GMX_HAVE_SIMD_UPDATE
269 /*! \brief Load (aligned) the contents of GMX_SIMD_REAL_WIDTH rvec elements sequentially into 3 SIMD registers
271 * The loaded output is:
272 * \p r0: { r[index][XX], r[index][YY], ... }
274 * \p r2: { ..., r[index+GMX_SIMD_REAL_WIDTH-1][YY], r[index+GMX_SIMD_REAL_WIDTH-1][ZZ] }
276 * \param[in] r Real to an rvec array, has to be aligned to SIMD register width
277 * \param[in] index Index of the first rvec triplet of reals to load
278 * \param[out] r0 Pointer to first SIMD register
279 * \param[out] r1 Pointer to second SIMD register
280 * \param[out] r2 Pointer to third SIMD register
282 static inline void simdLoadRvecs(const rvec
*r
,
288 const real
*realPtr
= r
[index
];
290 GMX_ASSERT(isSimdAligned(realPtr
), "Pointer should be SIMD aligned");
292 *r0
= simdLoad(realPtr
+ 0*GMX_SIMD_REAL_WIDTH
);
293 *r1
= simdLoad(realPtr
+ 1*GMX_SIMD_REAL_WIDTH
);
294 *r2
= simdLoad(realPtr
+ 2*GMX_SIMD_REAL_WIDTH
);
297 /*! \brief Store (aligned) 3 SIMD registers sequentially to GMX_SIMD_REAL_WIDTH rvec elements
299 * The stored output is:
300 * \p r[index] = { { r0[0], r0[1], ... }
302 * \p r[index+GMX_SIMD_REAL_WIDTH-1] = { ... , r2[GMX_SIMD_REAL_WIDTH-2], r2[GMX_SIMD_REAL_WIDTH-1] }
304 * \param[out] r Pointer to an rvec array, has to be aligned to SIMD register width
305 * \param[in] index Index of the first rvec triplet of reals to store to
306 * \param[in] r0 First SIMD register
307 * \param[in] r1 Second SIMD register
308 * \param[in] r2 Third SIMD register
310 static inline void simdStoreRvecs(rvec
*r
,
316 real
*realPtr
= r
[index
];
318 GMX_ASSERT(isSimdAligned(realPtr
), "Pointer should be SIMD aligned");
320 store(realPtr
+ 0*GMX_SIMD_REAL_WIDTH
, r0
);
321 store(realPtr
+ 1*GMX_SIMD_REAL_WIDTH
, r1
);
322 store(realPtr
+ 2*GMX_SIMD_REAL_WIDTH
, r2
);
325 /*! \brief Integrate using leap-frog with single group T-scaling and SIMD
327 * \param[in] start Index of first atom to update
328 * \param[in] nrend Last atom to update: \p nrend - 1
329 * \param[in] dt The time step
330 * \param[in] invMass 1/mass per atom
331 * \param[in] tcstat Temperature coupling information
332 * \param[in] x Input coordinates
333 * \param[out] xprime Updated coordinates
334 * \param[inout] v Velocities
335 * \param[in] f Forces
338 updateMDLeapfrogSimpleSimd(int start
,
341 const real
* gmx_restrict invMass
,
342 const t_grp_tcstat
* tcstat
,
343 const rvec
* gmx_restrict x
,
344 rvec
* gmx_restrict xprime
,
345 rvec
* gmx_restrict v
,
346 const rvec
* gmx_restrict f
)
348 SimdReal
timestep(dt
);
349 SimdReal
lambdaSystem(tcstat
[0].lambda
);
351 /* We declare variables here, since code is often slower when declaring them inside the loop */
353 /* Note: We should implement a proper PaddedVector, so we don't need this check */
354 GMX_ASSERT(isSimdAligned(invMass
), "invMass should be aligned");
356 for (int a
= start
; a
< nrend
; a
+= GMX_SIMD_REAL_WIDTH
)
358 SimdReal invMass0
, invMass1
, invMass2
;
359 expandScalarsToTriplets(simdLoad(invMass
+ a
),
360 &invMass0
, &invMass1
, &invMass2
);
364 simdLoadRvecs(v
, a
, &v0
, &v1
, &v2
);
365 simdLoadRvecs(f
, a
, &f0
, &f1
, &f2
);
367 v0
= fma(f0
*invMass0
, timestep
, lambdaSystem
*v0
);
368 v1
= fma(f1
*invMass1
, timestep
, lambdaSystem
*v1
);
369 v2
= fma(f2
*invMass2
, timestep
, lambdaSystem
*v2
);
371 simdStoreRvecs(v
, a
, v0
, v1
, v2
);
374 simdLoadRvecs(x
, a
, &x0
, &x1
, &x2
);
376 SimdReal xprime0
= fma(v0
, timestep
, x0
);
377 SimdReal xprime1
= fma(v1
, timestep
, x1
);
378 SimdReal xprime2
= fma(v2
, timestep
, x2
);
380 simdStoreRvecs(xprime
, a
, xprime0
, xprime1
, xprime2
);
384 #endif // GMX_HAVE_SIMD_UPDATE
386 /*! \brief Sets the NEMD acceleration type */
387 enum class AccelerationType
392 /*! \brief Integrate using leap-frog with support for everything
394 * \tparam accelerationType Type of NEMD acceleration
395 * \param[in] start Index of first atom to update
396 * \param[in] nrend Last atom to update: \p nrend - 1
397 * \param[in] doNoseHoover If to apply Nose-Hoover T-coupling
398 * \param[in] dt The time step
399 * \param[in] dtPressureCouple Time step for pressure coupling, is 0 when no pressure coupling should be applied at this step
400 * \param[in] ir The input parameter record
401 * \param[in] md Atom properties
402 * \param[in] ekind Kinetic energy data
403 * \param[in] box The box dimensions
404 * \param[in] x Input coordinates
405 * \param[out] xprime Updated coordinates
406 * \param[inout] v Velocities
407 * \param[in] f Forces
408 * \param[in] nh_vxi Nose-Hoover velocity scaling factors
409 * \param[in] M Parrinello-Rahman scaling matrix
411 template<AccelerationType accelerationType
>
413 updateMDLeapfrogGeneral(int start
,
417 real dtPressureCouple
,
418 const t_inputrec
* ir
,
419 const t_mdatoms
* md
,
420 const gmx_ekindata_t
* ekind
,
422 const rvec
* gmx_restrict x
,
423 rvec
* gmx_restrict xprime
,
424 rvec
* gmx_restrict v
,
425 const rvec
* gmx_restrict f
,
426 const double * gmx_restrict nh_vxi
,
429 /* This is a version of the leap-frog integrator that supports
430 * all combinations of T-coupling, P-coupling and NEMD.
431 * Nose-Hoover uses the reversible leap-frog integrator from
432 * Holian et al. Phys Rev E 52(3) : 2338, 1995
435 const unsigned short * cTC
= md
->cTC
;
436 const t_grp_tcstat
* tcstat
= ekind
->tcstat
;
438 const unsigned short * cACC
= md
->cACC
;
439 const rvec
* accel
= ir
->opts
.acc
;
440 const t_grp_acc
* grpstat
= ekind
->grpstat
;
442 const rvec
* gmx_restrict invMassPerDim
= md
->invMassPerDim
;
444 /* Initialize group values, changed later when multiple groups are used */
449 for (int n
= start
; n
< nrend
; n
++)
455 real lg
= tcstat
[gt
].lambda
;
458 real cosineZ
, vCosine
;
459 #ifdef __INTEL_COMPILER
460 #pragma warning( disable : 280 )
462 switch (accelerationType
)
464 case AccelerationType::none
:
465 copy_rvec(v
[n
], vRel
);
467 case AccelerationType::group
:
472 /* Avoid scaling the group velocity */
473 rvec_sub(v
[n
], grpstat
[ga
].u
, vRel
);
475 case AccelerationType::cosine
:
476 cosineZ
= std::cos(x
[n
][ZZ
]*static_cast<real
>(M_PI
)/box
[ZZ
][ZZ
]);
477 vCosine
= cosineZ
*ekind
->cosacc
.vcos
;
478 /* Avoid scaling the cosine profile velocity */
479 copy_rvec(v
[n
], vRel
);
486 /* Here we account for multiple time stepping, by increasing
487 * the Nose-Hoover correction by nsttcouple
489 factorNH
= 0.5*ir
->nsttcouple
*dt
*nh_vxi
[gt
];
492 for (int d
= 0; d
< DIM
; d
++)
495 (lg
*vRel
[d
] + (f
[n
][d
]*invMassPerDim
[n
][d
]*dt
- factorNH
*vRel
[d
]
496 - dtPressureCouple
*iprod(M
[d
], vRel
)))/(1 + factorNH
);
497 switch (accelerationType
)
499 case AccelerationType::none
:
501 case AccelerationType::group
:
502 /* Add back the mean velocity and apply acceleration */
503 vNew
+= grpstat
[ga
].u
[d
] + accel
[ga
][d
]*dt
;
505 case AccelerationType::cosine
:
508 /* Add back the mean velocity and apply acceleration */
509 vNew
+= vCosine
+ cosineZ
*ekind
->cosacc
.cos_accel
*dt
;
514 xprime
[n
][d
] = x
[n
][d
] + vNew
*dt
;
519 /*! \brief Handles the Leap-frog MD x and v integration */
520 static void do_update_md(int start
,
524 const t_inputrec
* ir
,
525 const t_mdatoms
* md
,
526 const gmx_ekindata_t
* ekind
,
528 const rvec
* gmx_restrict x
,
529 rvec
* gmx_restrict xprime
,
530 rvec
* gmx_restrict v
,
531 const rvec
* gmx_restrict f
,
532 const double * gmx_restrict nh_vxi
,
535 GMX_ASSERT(nrend
== start
|| xprime
!= x
, "For SIMD optimization certain compilers need to have xprime != x");
537 /* Note: Berendsen pressure scaling is handled after do_update_md() */
538 bool doTempCouple
= isTemperatureCouplingStep(step
, ir
);
539 bool doNoseHoover
= (ir
->etc
== etcNOSEHOOVER
&& doTempCouple
);
540 bool doParrinelloRahman
= (ir
->epc
== epcPARRINELLORAHMAN
&& isPressureCouplingStep(step
, ir
));
541 bool doPROffDiagonal
= (doParrinelloRahman
&& (M
[YY
][XX
] != 0 || M
[ZZ
][XX
] != 0 || M
[ZZ
][YY
] != 0));
543 real dtPressureCouple
= (doParrinelloRahman
? ir
->nstpcouple
*dt
: 0);
545 /* NEMD (also cosine) acceleration is applied in updateMDLeapFrogGeneral */
546 bool doAcceleration
= (ekind
->bNEMD
|| ekind
->cosacc
.cos_accel
!= 0);
548 if (doNoseHoover
|| doPROffDiagonal
|| doAcceleration
)
551 if (!doParrinelloRahman
)
553 /* We should not apply PR scaling at this step */
563 updateMDLeapfrogGeneral
<AccelerationType::none
>
564 (start
, nrend
, doNoseHoover
, dt
, dtPressureCouple
,
565 ir
, md
, ekind
, box
, x
, xprime
, v
, f
, nh_vxi
, stepM
);
567 else if (ekind
->bNEMD
)
569 updateMDLeapfrogGeneral
<AccelerationType::group
>
570 (start
, nrend
, doNoseHoover
, dt
, dtPressureCouple
,
571 ir
, md
, ekind
, box
, x
, xprime
, v
, f
, nh_vxi
, stepM
);
575 updateMDLeapfrogGeneral
<AccelerationType::cosine
>
576 (start
, nrend
, doNoseHoover
, dt
, dtPressureCouple
,
577 ir
, md
, ekind
, box
, x
, xprime
, v
, f
, nh_vxi
, stepM
);
582 /* Use a simple and thus more efficient integration loop. */
583 /* The simple loop does not check for particle type (so it can
584 * be vectorized), which means we need to clear the velocities
585 * of virtual sites in advance, when present. Note that vsite
586 * velocities are computed after update and constraints from
587 * their displacement.
591 /* Note: The overhead of this loop is completely neligible */
592 clearVsiteVelocities(start
, nrend
, md
->ptype
, v
);
595 /* We determine if we have a single T-coupling lambda value for all
596 * atoms. That allows for better SIMD acceleration in the template.
597 * If we do not do temperature coupling (in the run or this step),
598 * all scaling values are 1, so we effectively have a single value.
600 bool haveSingleTempScaleValue
= (!doTempCouple
|| ekind
->ngtc
== 1);
602 /* Extract some pointers needed by all cases */
603 const unsigned short *cTC
= md
->cTC
;
604 const t_grp_tcstat
*tcstat
= ekind
->tcstat
;
605 const rvec
*invMassPerDim
= md
->invMassPerDim
;
607 if (doParrinelloRahman
)
609 GMX_ASSERT(!doPROffDiagonal
, "updateMDLeapfrogSimple only support diagonal Parrinello-Rahman scaling matrices");
612 for (int d
= 0; d
< DIM
; d
++)
617 if (haveSingleTempScaleValue
)
619 updateMDLeapfrogSimple
620 <NumTempScaleValues::single
,
621 ApplyParrinelloRahmanVScaling::diagonal
>
622 (start
, nrend
, dt
, dtPressureCouple
,
623 invMassPerDim
, tcstat
, cTC
, diagM
, x
, xprime
, v
, f
);
627 updateMDLeapfrogSimple
628 <NumTempScaleValues::multiple
,
629 ApplyParrinelloRahmanVScaling::diagonal
>
630 (start
, nrend
, dt
, dtPressureCouple
,
631 invMassPerDim
, tcstat
, cTC
, diagM
, x
, xprime
, v
, f
);
636 if (haveSingleTempScaleValue
)
638 /* Note that modern compilers are pretty good at vectorizing
639 * updateMDLeapfrogSimple(). But the SIMD version will still
640 * be faster because invMass lowers the cache pressure
641 * compared to invMassPerDim.
643 #if GMX_HAVE_SIMD_UPDATE
644 /* Check if we can use invmass instead of invMassPerDim */
645 if (!md
->havePartiallyFrozenAtoms
)
647 updateMDLeapfrogSimpleSimd
648 (start
, nrend
, dt
, md
->invmass
, tcstat
, x
, xprime
, v
, f
);
653 updateMDLeapfrogSimple
654 <NumTempScaleValues::single
,
655 ApplyParrinelloRahmanVScaling::no
>
656 (start
, nrend
, dt
, dtPressureCouple
,
657 invMassPerDim
, tcstat
, cTC
, nullptr, x
, xprime
, v
, f
);
662 updateMDLeapfrogSimple
663 <NumTempScaleValues::multiple
,
664 ApplyParrinelloRahmanVScaling::no
>
665 (start
, nrend
, dt
, dtPressureCouple
,
666 invMassPerDim
, tcstat
, cTC
, nullptr, x
, xprime
, v
, f
);
672 static void do_update_vv_vel(int start
, int nrend
, real dt
,
673 rvec accel
[], ivec nFreeze
[], real invmass
[],
674 unsigned short ptype
[], unsigned short cFREEZE
[],
675 unsigned short cACC
[], rvec v
[], const rvec f
[],
676 gmx_bool bExtended
, real veta
, real alpha
)
684 g
= 0.25*dt
*veta
*alpha
;
686 mv2
= gmx::series_sinhx(g
);
693 for (n
= start
; n
< nrend
; n
++)
695 real w_dt
= invmass
[n
]*dt
;
705 for (d
= 0; d
< DIM
; d
++)
707 if ((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
) && !nFreeze
[gf
][d
])
709 v
[n
][d
] = mv1
*(mv1
*v
[n
][d
] + 0.5*(w_dt
*mv2
*f
[n
][d
]))+0.5*accel
[ga
][d
]*dt
;
717 } /* do_update_vv_vel */
719 static void do_update_vv_pos(int start
, int nrend
, real dt
,
721 unsigned short ptype
[], unsigned short cFREEZE
[],
722 const rvec x
[], rvec xprime
[], rvec v
[],
723 gmx_bool bExtended
, real veta
)
729 /* Would it make more sense if Parrinello-Rahman was put here? */
734 mr2
= gmx::series_sinhx(g
);
742 for (n
= start
; n
< nrend
; n
++)
750 for (d
= 0; d
< DIM
; d
++)
752 if ((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
) && !nFreeze
[gf
][d
])
754 xprime
[n
][d
] = mr1
*(mr1
*x
[n
][d
]+mr2
*dt
*v
[n
][d
]);
758 xprime
[n
][d
] = x
[n
][d
];
762 } /* do_update_vv_pos */
764 static gmx_stochd_t
*init_stochd(const t_inputrec
*ir
)
770 const t_grpopts
*opts
= &ir
->opts
;
771 int ngtc
= opts
->ngtc
;
775 snew(sd
->bd_rf
, ngtc
);
777 else if (EI_SD(ir
->eI
))
780 snew(sd
->sdsig
, ngtc
);
782 gmx_sd_const_t
*sdc
= sd
->sdc
;
784 for (int gt
= 0; gt
< ngtc
; gt
++)
786 if (opts
->tau_t
[gt
] > 0)
788 sdc
[gt
].em
= std::exp(-ir
->delta_t
/opts
->tau_t
[gt
]);
792 /* No friction and noise on this group */
797 else if (ETC_ANDERSEN(ir
->etc
))
799 snew(sd
->randomize_group
, ngtc
);
800 snew(sd
->boltzfac
, ngtc
);
802 /* for now, assume that all groups, if randomized, are randomized at the same rate, i.e. tau_t is the same. */
803 /* since constraint groups don't necessarily match up with temperature groups! This is checked in readir.c */
805 for (int gt
= 0; gt
< ngtc
; gt
++)
807 real reft
= std::max
<real
>(0, opts
->ref_t
[gt
]);
808 if ((opts
->tau_t
[gt
] > 0) && (reft
> 0)) /* tau_t or ref_t = 0 means that no randomization is done */
810 sd
->randomize_group
[gt
] = TRUE
;
811 sd
->boltzfac
[gt
] = BOLTZ
*opts
->ref_t
[gt
];
815 sd
->randomize_group
[gt
] = FALSE
;
823 void update_temperature_constants(gmx_update_t
*upd
, const t_inputrec
*ir
)
827 if (ir
->bd_fric
!= 0)
829 for (int gt
= 0; gt
< ir
->opts
.ngtc
; gt
++)
831 upd
->sd
->bd_rf
[gt
] = std::sqrt(2.0*BOLTZ
*ir
->opts
.ref_t
[gt
]/(ir
->bd_fric
*ir
->delta_t
));
836 for (int gt
= 0; gt
< ir
->opts
.ngtc
; gt
++)
838 upd
->sd
->bd_rf
[gt
] = std::sqrt(2.0*BOLTZ
*ir
->opts
.ref_t
[gt
]);
844 for (int gt
= 0; gt
< ir
->opts
.ngtc
; gt
++)
846 real kT
= BOLTZ
*ir
->opts
.ref_t
[gt
];
847 /* The mass is accounted for later, since this differs per atom */
848 upd
->sd
->sdsig
[gt
].V
= std::sqrt(kT
*(1 - upd
->sd
->sdc
[gt
].em
*upd
->sd
->sdc
[gt
].em
));
853 gmx_update_t
*init_update(const t_inputrec
*ir
)
855 gmx_update_t
*upd
= new(gmx_update_t
);
857 if (ir
->eI
== eiBD
|| EI_SD(ir
->eI
) || ir
->etc
== etcVRESCALE
|| ETC_ANDERSEN(ir
->etc
))
859 upd
->sd
= init_stochd(ir
);
862 update_temperature_constants(upd
, ir
);
869 void update_realloc(gmx_update_t
*upd
, int natoms
)
871 GMX_ASSERT(upd
, "upd must be allocated before its fields can be reallocated");
873 upd
->xp
.resize(gmx::paddedRVecVectorSize(natoms
));
876 /*! \brief Sets the SD update type */
877 enum class SDUpdate
: int
879 ForcesOnly
, FrictionAndNoiseOnly
, Combined
882 /*! \brief SD integrator update
884 * Two phases are required in the general case of a constrained
885 * update, the first phase from the contribution of forces, before
886 * applying constraints, and then a second phase applying the friction
887 * and noise, and then further constraining. For details, see
890 * Without constraints, the two phases can be combined, for
893 * Thus three instantiations of this templated function will be made,
894 * two with only one contribution, and one with both contributions. */
895 template <SDUpdate updateType
>
897 doSDUpdateGeneral(gmx_stochd_t
*sd
,
898 int start
, int nrend
, real dt
,
899 rvec accel
[], ivec nFreeze
[],
900 real invmass
[], unsigned short ptype
[],
901 unsigned short cFREEZE
[], unsigned short cACC
[],
902 unsigned short cTC
[],
903 const rvec x
[], rvec xprime
[], rvec v
[], const rvec f
[],
904 gmx_int64_t step
, int seed
, int *gatindex
)
906 if (updateType
!= SDUpdate::FrictionAndNoiseOnly
)
908 GMX_ASSERT(f
!= nullptr, "SD update with forces requires forces");
909 GMX_ASSERT(cACC
!= nullptr, "SD update with forces requires acceleration groups");
911 if (updateType
!= SDUpdate::ForcesOnly
)
913 GMX_ASSERT(cTC
!= nullptr, "SD update with noise requires temperature groups");
916 // Even 0 bits internal counter gives 2x64 ints (more than enough for three table lookups)
917 gmx::ThreeFry2x64
<0> rng(seed
, gmx::RandomDomain::UpdateCoordinates
);
918 gmx::TabulatedNormalDistribution
<real
, 14> dist
;
920 gmx_sd_const_t
*sdc
= sd
->sdc
;
921 gmx_sd_sigma_t
*sig
= sd
->sdsig
;
923 for (int n
= start
; n
< nrend
; n
++)
925 int globalAtomIndex
= gatindex
? gatindex
[n
] : n
;
926 rng
.restart(step
, globalAtomIndex
);
929 real inverseMass
= invmass
[n
];
930 real invsqrtMass
= std::sqrt(inverseMass
);
932 int freezeGroup
= cFREEZE
? cFREEZE
[n
] : 0;
933 int accelerationGroup
= cACC
? cACC
[n
] : 0;
934 int temperatureGroup
= cTC
? cTC
[n
] : 0;
936 for (int d
= 0; d
< DIM
; d
++)
938 if ((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
) && !nFreeze
[freezeGroup
][d
])
940 if (updateType
== SDUpdate::ForcesOnly
)
942 real vn
= v
[n
][d
] + (inverseMass
*f
[n
][d
] + accel
[accelerationGroup
][d
])*dt
;
944 // Simple position update.
945 xprime
[n
][d
] = x
[n
][d
] + v
[n
][d
]*dt
;
947 else if (updateType
== SDUpdate::FrictionAndNoiseOnly
)
950 v
[n
][d
] = (vn
*sdc
[temperatureGroup
].em
+
951 invsqrtMass
*sig
[temperatureGroup
].V
*dist(rng
));
952 // The previous phase already updated the
953 // positions with a full v*dt term that must
954 // now be half removed.
955 xprime
[n
][d
] = xprime
[n
][d
] + 0.5*(v
[n
][d
] - vn
)*dt
;
959 real vn
= v
[n
][d
] + (inverseMass
*f
[n
][d
] + accel
[accelerationGroup
][d
])*dt
;
960 v
[n
][d
] = (vn
*sdc
[temperatureGroup
].em
+
961 invsqrtMass
*sig
[temperatureGroup
].V
*dist(rng
));
962 // Here we include half of the friction+noise
963 // update of v into the position update.
964 xprime
[n
][d
] = x
[n
][d
] + 0.5*(vn
+ v
[n
][d
])*dt
;
969 // When using constraints, the update is split into
970 // two phases, but we only need to zero the update of
971 // virtual, shell or frozen particles in at most one
973 if (updateType
!= SDUpdate::FrictionAndNoiseOnly
)
976 xprime
[n
][d
] = x
[n
][d
];
983 static void do_update_bd(int start
, int nrend
, real dt
,
985 real invmass
[], unsigned short ptype
[],
986 unsigned short cFREEZE
[], unsigned short cTC
[],
987 const rvec x
[], rvec xprime
[], rvec v
[],
988 const rvec f
[], real friction_coefficient
,
989 real
*rf
, gmx_int64_t step
, int seed
,
992 /* note -- these appear to be full step velocities . . . */
997 // Use 1 bit of internal counters to give us 2*2 64-bits values per stream
998 // Each 64-bit value is enough for 4 normal distribution table numbers.
999 gmx::ThreeFry2x64
<0> rng(seed
, gmx::RandomDomain::UpdateCoordinates
);
1000 gmx::TabulatedNormalDistribution
<real
, 14> dist
;
1002 if (friction_coefficient
!= 0)
1004 invfr
= 1.0/friction_coefficient
;
1007 for (n
= start
; (n
< nrend
); n
++)
1009 int ng
= gatindex
? gatindex
[n
] : n
;
1011 rng
.restart(step
, ng
);
1022 for (d
= 0; (d
< DIM
); d
++)
1024 if ((ptype
[n
] != eptVSite
) && (ptype
[n
] != eptShell
) && !nFreeze
[gf
][d
])
1026 if (friction_coefficient
!= 0)
1028 vn
= invfr
*f
[n
][d
] + rf
[gt
]*dist(rng
);
1032 /* NOTE: invmass = 2/(mass*friction_constant*dt) */
1033 vn
= 0.5*invmass
[n
]*f
[n
][d
]*dt
1034 + std::sqrt(0.5*invmass
[n
])*rf
[gt
]*dist(rng
);
1038 xprime
[n
][d
] = x
[n
][d
]+vn
*dt
;
1043 xprime
[n
][d
] = x
[n
][d
];
1049 static void calc_ke_part_normal(rvec v
[], t_grpopts
*opts
, t_mdatoms
*md
,
1050 gmx_ekindata_t
*ekind
, t_nrnb
*nrnb
, gmx_bool bEkinAveVel
)
1053 t_grp_tcstat
*tcstat
= ekind
->tcstat
;
1054 t_grp_acc
*grpstat
= ekind
->grpstat
;
1055 int nthread
, thread
;
1057 /* three main: VV with AveVel, vv with AveEkin, leap with AveEkin. Leap with AveVel is also
1058 an option, but not supported now.
1059 bEkinAveVel: If TRUE, we sum into ekin, if FALSE, into ekinh.
1062 /* group velocities are calculated in update_ekindata and
1063 * accumulated in acumulate_groups.
1064 * Now the partial global and groups ekin.
1066 for (g
= 0; (g
< opts
->ngtc
); g
++)
1068 copy_mat(tcstat
[g
].ekinh
, tcstat
[g
].ekinh_old
);
1071 clear_mat(tcstat
[g
].ekinf
);
1072 tcstat
[g
].ekinscalef_nhc
= 1.0; /* need to clear this -- logic is complicated! */
1076 clear_mat(tcstat
[g
].ekinh
);
1079 ekind
->dekindl_old
= ekind
->dekindl
;
1080 nthread
= gmx_omp_nthreads_get(emntUpdate
);
1082 #pragma omp parallel for num_threads(nthread) schedule(static)
1083 for (thread
= 0; thread
< nthread
; thread
++)
1085 // This OpenMP only loops over arrays and does not call any functions
1086 // or memory allocation. It should not be able to throw, so for now
1087 // we do not need a try/catch wrapper.
1088 int start_t
, end_t
, n
;
1096 start_t
= ((thread
+0)*md
->homenr
)/nthread
;
1097 end_t
= ((thread
+1)*md
->homenr
)/nthread
;
1099 ekin_sum
= ekind
->ekin_work
[thread
];
1100 dekindl_sum
= ekind
->dekindl_work
[thread
];
1102 for (gt
= 0; gt
< opts
->ngtc
; gt
++)
1104 clear_mat(ekin_sum
[gt
]);
1110 for (n
= start_t
; n
< end_t
; n
++)
1120 hm
= 0.5*md
->massT
[n
];
1122 for (d
= 0; (d
< DIM
); d
++)
1124 v_corrt
[d
] = v
[n
][d
] - grpstat
[ga
].u
[d
];
1126 for (d
= 0; (d
< DIM
); d
++)
1128 for (m
= 0; (m
< DIM
); m
++)
1130 /* if we're computing a full step velocity, v_corrt[d] has v(t). Otherwise, v(t+dt/2) */
1131 ekin_sum
[gt
][m
][d
] += hm
*v_corrt
[m
]*v_corrt
[d
];
1134 if (md
->nMassPerturbed
&& md
->bPerturbed
[n
])
1137 0.5*(md
->massB
[n
] - md
->massA
[n
])*iprod(v_corrt
, v_corrt
);
1143 for (thread
= 0; thread
< nthread
; thread
++)
1145 for (g
= 0; g
< opts
->ngtc
; g
++)
1149 m_add(tcstat
[g
].ekinf
, ekind
->ekin_work
[thread
][g
],
1154 m_add(tcstat
[g
].ekinh
, ekind
->ekin_work
[thread
][g
],
1159 ekind
->dekindl
+= *ekind
->dekindl_work
[thread
];
1162 inc_nrnb(nrnb
, eNR_EKIN
, md
->homenr
);
1165 static void calc_ke_part_visc(matrix box
, rvec x
[], rvec v
[],
1166 t_grpopts
*opts
, t_mdatoms
*md
,
1167 gmx_ekindata_t
*ekind
,
1168 t_nrnb
*nrnb
, gmx_bool bEkinAveVel
)
1170 int start
= 0, homenr
= md
->homenr
;
1171 int g
, d
, n
, m
, gt
= 0;
1174 t_grp_tcstat
*tcstat
= ekind
->tcstat
;
1175 t_cos_acc
*cosacc
= &(ekind
->cosacc
);
1180 for (g
= 0; g
< opts
->ngtc
; g
++)
1182 copy_mat(ekind
->tcstat
[g
].ekinh
, ekind
->tcstat
[g
].ekinh_old
);
1183 clear_mat(ekind
->tcstat
[g
].ekinh
);
1185 ekind
->dekindl_old
= ekind
->dekindl
;
1187 fac
= 2*M_PI
/box
[ZZ
][ZZ
];
1190 for (n
= start
; n
< start
+homenr
; n
++)
1196 hm
= 0.5*md
->massT
[n
];
1198 /* Note that the times of x and v differ by half a step */
1199 /* MRS -- would have to be changed for VV */
1200 cosz
= std::cos(fac
*x
[n
][ZZ
]);
1201 /* Calculate the amplitude of the new velocity profile */
1202 mvcos
+= 2*cosz
*md
->massT
[n
]*v
[n
][XX
];
1204 copy_rvec(v
[n
], v_corrt
);
1205 /* Subtract the profile for the kinetic energy */
1206 v_corrt
[XX
] -= cosz
*cosacc
->vcos
;
1207 for (d
= 0; (d
< DIM
); d
++)
1209 for (m
= 0; (m
< DIM
); m
++)
1211 /* if we're computing a full step velocity, v_corrt[d] has v(t). Otherwise, v(t+dt/2) */
1214 tcstat
[gt
].ekinf
[m
][d
] += hm
*v_corrt
[m
]*v_corrt
[d
];
1218 tcstat
[gt
].ekinh
[m
][d
] += hm
*v_corrt
[m
]*v_corrt
[d
];
1222 if (md
->nPerturbed
&& md
->bPerturbed
[n
])
1224 /* The minus sign here might be confusing.
1225 * The kinetic contribution from dH/dl doesn't come from
1226 * d m(l)/2 v^2 / dl, but rather from d p^2/2m(l) / dl,
1227 * where p are the momenta. The difference is only a minus sign.
1229 dekindl
-= 0.5*(md
->massB
[n
] - md
->massA
[n
])*iprod(v_corrt
, v_corrt
);
1232 ekind
->dekindl
= dekindl
;
1233 cosacc
->mvcos
= mvcos
;
1235 inc_nrnb(nrnb
, eNR_EKIN
, homenr
);
1238 void calc_ke_part(t_state
*state
, t_grpopts
*opts
, t_mdatoms
*md
,
1239 gmx_ekindata_t
*ekind
, t_nrnb
*nrnb
, gmx_bool bEkinAveVel
)
1241 if (ekind
->cosacc
.cos_accel
== 0)
1243 calc_ke_part_normal(as_rvec_array(state
->v
.data()), opts
, md
, ekind
, nrnb
, bEkinAveVel
);
1247 calc_ke_part_visc(state
->box
, as_rvec_array(state
->x
.data()), as_rvec_array(state
->v
.data()), opts
, md
, ekind
, nrnb
, bEkinAveVel
);
1251 extern void init_ekinstate(ekinstate_t
*ekinstate
, const t_inputrec
*ir
)
1253 ekinstate
->ekin_n
= ir
->opts
.ngtc
;
1254 snew(ekinstate
->ekinh
, ekinstate
->ekin_n
);
1255 snew(ekinstate
->ekinf
, ekinstate
->ekin_n
);
1256 snew(ekinstate
->ekinh_old
, ekinstate
->ekin_n
);
1257 ekinstate
->ekinscalef_nhc
.resize(ekinstate
->ekin_n
);
1258 ekinstate
->ekinscaleh_nhc
.resize(ekinstate
->ekin_n
);
1259 ekinstate
->vscale_nhc
.resize(ekinstate
->ekin_n
);
1260 ekinstate
->dekindl
= 0;
1261 ekinstate
->mvcos
= 0;
1264 void update_ekinstate(ekinstate_t
*ekinstate
, gmx_ekindata_t
*ekind
)
1268 for (i
= 0; i
< ekinstate
->ekin_n
; i
++)
1270 copy_mat(ekind
->tcstat
[i
].ekinh
, ekinstate
->ekinh
[i
]);
1271 copy_mat(ekind
->tcstat
[i
].ekinf
, ekinstate
->ekinf
[i
]);
1272 copy_mat(ekind
->tcstat
[i
].ekinh_old
, ekinstate
->ekinh_old
[i
]);
1273 ekinstate
->ekinscalef_nhc
[i
] = ekind
->tcstat
[i
].ekinscalef_nhc
;
1274 ekinstate
->ekinscaleh_nhc
[i
] = ekind
->tcstat
[i
].ekinscaleh_nhc
;
1275 ekinstate
->vscale_nhc
[i
] = ekind
->tcstat
[i
].vscale_nhc
;
1278 copy_mat(ekind
->ekin
, ekinstate
->ekin_total
);
1279 ekinstate
->dekindl
= ekind
->dekindl
;
1280 ekinstate
->mvcos
= ekind
->cosacc
.mvcos
;
1284 void restore_ekinstate_from_state(const t_commrec
*cr
,
1285 gmx_ekindata_t
*ekind
, const ekinstate_t
*ekinstate
)
1291 for (i
= 0; i
< ekinstate
->ekin_n
; i
++)
1293 copy_mat(ekinstate
->ekinh
[i
], ekind
->tcstat
[i
].ekinh
);
1294 copy_mat(ekinstate
->ekinf
[i
], ekind
->tcstat
[i
].ekinf
);
1295 copy_mat(ekinstate
->ekinh_old
[i
], ekind
->tcstat
[i
].ekinh_old
);
1296 ekind
->tcstat
[i
].ekinscalef_nhc
= ekinstate
->ekinscalef_nhc
[i
];
1297 ekind
->tcstat
[i
].ekinscaleh_nhc
= ekinstate
->ekinscaleh_nhc
[i
];
1298 ekind
->tcstat
[i
].vscale_nhc
= ekinstate
->vscale_nhc
[i
];
1301 copy_mat(ekinstate
->ekin_total
, ekind
->ekin
);
1303 ekind
->dekindl
= ekinstate
->dekindl
;
1304 ekind
->cosacc
.mvcos
= ekinstate
->mvcos
;
1305 n
= ekinstate
->ekin_n
;
1310 gmx_bcast(sizeof(n
), &n
, cr
);
1311 for (i
= 0; i
< n
; i
++)
1313 gmx_bcast(DIM
*DIM
*sizeof(ekind
->tcstat
[i
].ekinh
[0][0]),
1314 ekind
->tcstat
[i
].ekinh
[0], cr
);
1315 gmx_bcast(DIM
*DIM
*sizeof(ekind
->tcstat
[i
].ekinf
[0][0]),
1316 ekind
->tcstat
[i
].ekinf
[0], cr
);
1317 gmx_bcast(DIM
*DIM
*sizeof(ekind
->tcstat
[i
].ekinh_old
[0][0]),
1318 ekind
->tcstat
[i
].ekinh_old
[0], cr
);
1320 gmx_bcast(sizeof(ekind
->tcstat
[i
].ekinscalef_nhc
),
1321 &(ekind
->tcstat
[i
].ekinscalef_nhc
), cr
);
1322 gmx_bcast(sizeof(ekind
->tcstat
[i
].ekinscaleh_nhc
),
1323 &(ekind
->tcstat
[i
].ekinscaleh_nhc
), cr
);
1324 gmx_bcast(sizeof(ekind
->tcstat
[i
].vscale_nhc
),
1325 &(ekind
->tcstat
[i
].vscale_nhc
), cr
);
1327 gmx_bcast(DIM
*DIM
*sizeof(ekind
->ekin
[0][0]),
1328 ekind
->ekin
[0], cr
);
1330 gmx_bcast(sizeof(ekind
->dekindl
), &ekind
->dekindl
, cr
);
1331 gmx_bcast(sizeof(ekind
->cosacc
.mvcos
), &ekind
->cosacc
.mvcos
, cr
);
1335 void set_deform_reference_box(gmx_update_t
*upd
, gmx_int64_t step
, matrix box
)
1337 upd
->deformref_step
= step
;
1338 copy_mat(box
, upd
->deformref_box
);
1341 static void deform(gmx_update_t
*upd
,
1342 int start
, int homenr
, rvec x
[], matrix box
,
1343 const t_inputrec
*ir
, gmx_int64_t step
)
1345 matrix bnew
, invbox
, mu
;
1349 elapsed_time
= (step
+ 1 - upd
->deformref_step
)*ir
->delta_t
;
1350 copy_mat(box
, bnew
);
1351 for (i
= 0; i
< DIM
; i
++)
1353 for (j
= 0; j
< DIM
; j
++)
1355 if (ir
->deform
[i
][j
] != 0)
1358 upd
->deformref_box
[i
][j
] + elapsed_time
*ir
->deform
[i
][j
];
1362 /* We correct the off-diagonal elements,
1363 * which can grow indefinitely during shearing,
1364 * so the shifts do not get messed up.
1366 for (i
= 1; i
< DIM
; i
++)
1368 for (j
= i
-1; j
>= 0; j
--)
1370 while (bnew
[i
][j
] - box
[i
][j
] > 0.5*bnew
[j
][j
])
1372 rvec_dec(bnew
[i
], bnew
[j
]);
1374 while (bnew
[i
][j
] - box
[i
][j
] < -0.5*bnew
[j
][j
])
1376 rvec_inc(bnew
[i
], bnew
[j
]);
1380 gmx::invertBoxMatrix(box
, invbox
);
1381 copy_mat(bnew
, box
);
1382 mmul_ur0(box
, invbox
, mu
);
1384 for (i
= start
; i
< start
+homenr
; i
++)
1386 x
[i
][XX
] = mu
[XX
][XX
]*x
[i
][XX
]+mu
[YY
][XX
]*x
[i
][YY
]+mu
[ZZ
][XX
]*x
[i
][ZZ
];
1387 x
[i
][YY
] = mu
[YY
][YY
]*x
[i
][YY
]+mu
[ZZ
][YY
]*x
[i
][ZZ
];
1388 x
[i
][ZZ
] = mu
[ZZ
][ZZ
]*x
[i
][ZZ
];
1392 void update_tcouple(gmx_int64_t step
,
1393 t_inputrec
*inputrec
,
1395 gmx_ekindata_t
*ekind
,
1400 bool doTemperatureCoupling
= false;
1402 /* if using vv with trotter decomposition methods, we do this elsewhere in the code */
1403 if (!(EI_VV(inputrec
->eI
) &&
1404 (inputrecNvtTrotter(inputrec
) || inputrecNptTrotter(inputrec
) || inputrecNphTrotter(inputrec
))))
1406 doTemperatureCoupling
= isTemperatureCouplingStep(step
, inputrec
);
1409 if (doTemperatureCoupling
)
1411 real dttc
= inputrec
->nsttcouple
*inputrec
->delta_t
;
1413 switch (inputrec
->etc
)
1418 berendsen_tcoupl(inputrec
, ekind
, dttc
, state
->therm_integral
);
1421 nosehoover_tcoupl(&(inputrec
->opts
), ekind
, dttc
,
1422 state
->nosehoover_xi
.data(), state
->nosehoover_vxi
.data(), MassQ
);
1425 vrescale_tcoupl(inputrec
, step
, ekind
, dttc
,
1426 state
->therm_integral
.data());
1429 /* rescale in place here */
1430 if (EI_VV(inputrec
->eI
))
1432 rescale_velocities(ekind
, md
, 0, md
->homenr
, as_rvec_array(state
->v
.data()));
1437 /* Set the T scaling lambda to 1 to have no scaling */
1438 for (int i
= 0; (i
< inputrec
->opts
.ngtc
); i
++)
1440 ekind
->tcstat
[i
].lambda
= 1.0;
1445 /*! \brief Computes the atom range for a thread to operate on, ensured SIMD aligned ranges
1447 * \param[in] numThreads The number of threads to divide atoms over
1448 * \param[in] threadIndex The thread to get the range for
1449 * \param[in] numAtoms The total number of atoms (on this rank)
1450 * \param[out] startAtom The start of the atom range
1451 * \param[out] endAtom The end of the atom range, note that this is in general not a multiple of the SIMD width
1453 static void getThreadAtomRange(int numThreads
, int threadIndex
, int numAtoms
,
1454 int *startAtom
, int *endAtom
)
1456 #if GMX_HAVE_SIMD_UPDATE
1457 constexpr int blockSize
= GMX_SIMD_REAL_WIDTH
;
1459 constexpr int blockSize
= 1;
1461 int numBlocks
= (numAtoms
+ blockSize
- 1)/blockSize
;
1463 *startAtom
= ((numBlocks
* threadIndex
)/numThreads
)*blockSize
;
1464 *endAtom
= ((numBlocks
*(threadIndex
+ 1))/numThreads
)*blockSize
;
1465 if (threadIndex
== numThreads
- 1)
1467 *endAtom
= numAtoms
;
1471 void update_pcouple_before_coordinates(FILE *fplog
,
1473 const t_inputrec
*inputrec
,
1475 matrix parrinellorahmanMu
,
1479 /* Berendsen P-coupling is completely handled after the coordinate update.
1480 * Trotter P-coupling is handled by separate calls to trotter_update().
1482 if (inputrec
->epc
== epcPARRINELLORAHMAN
&&
1483 isPressureCouplingStep(step
, inputrec
))
1485 real dtpc
= inputrec
->nstpcouple
*inputrec
->delta_t
;
1487 parrinellorahman_pcoupl(fplog
, step
, inputrec
, dtpc
, state
->pres_prev
,
1488 state
->box
, state
->box_rel
, state
->boxv
,
1489 M
, parrinellorahmanMu
, bInitStep
);
1493 void constrain_velocities(gmx_int64_t step
,
1494 real
*dvdlambda
, /* the contribution to be added to the bonded interactions */
1495 const t_inputrec
*inputrec
, /* input record and box stuff */
1501 const t_commrec
*cr
,
1502 const gmx_multisim_t
*ms
,
1504 gmx_wallcycle_t wcycle
,
1505 gmx::Constraints
*constr
,
1517 * APPLY CONSTRAINTS:
1524 /* clear out constraints before applying */
1525 clear_mat(vir_part
);
1527 /* Constrain the coordinates upd->xp */
1528 wallcycle_start(wcycle
, ewcCONSTR
);
1530 constrain(nullptr, do_log
, do_ene
, constr
, idef
,
1531 inputrec
, cr
, ms
, step
, 1, 1.0, md
,
1532 as_rvec_array(state
->x
.data()), as_rvec_array(state
->v
.data()), as_rvec_array(state
->v
.data()),
1533 bMolPBC
, state
->box
,
1534 state
->lambda
[efptBONDED
], dvdlambda
,
1535 nullptr, bCalcVir
? &vir_con
: nullptr, nrnb
, econqVeloc
);
1537 wallcycle_stop(wcycle
, ewcCONSTR
);
1541 m_add(vir_part
, vir_con
, vir_part
);
1546 void constrain_coordinates(gmx_int64_t step
,
1547 real
*dvdlambda
, /* the contribution to be added to the bonded interactions */
1548 const t_inputrec
*inputrec
, /* input record and box stuff */
1554 const t_commrec
*cr
,
1555 const gmx_multisim_t
*ms
,
1557 gmx_wallcycle_t wcycle
,
1559 gmx::Constraints
*constr
,
1572 /* clear out constraints before applying */
1573 clear_mat(vir_part
);
1575 /* Constrain the coordinates upd->xp */
1576 wallcycle_start(wcycle
, ewcCONSTR
);
1578 constrain(nullptr, do_log
, do_ene
, constr
, idef
,
1579 inputrec
, cr
, ms
, step
, 1, 1.0, md
,
1580 as_rvec_array(state
->x
.data()), as_rvec_array(upd
->xp
.data()), nullptr,
1581 bMolPBC
, state
->box
,
1582 state
->lambda
[efptBONDED
], dvdlambda
,
1583 as_rvec_array(state
->v
.data()), bCalcVir
? &vir_con
: nullptr, nrnb
, econqCoord
);
1585 wallcycle_stop(wcycle
, ewcCONSTR
);
1589 m_add(vir_part
, vir_con
, vir_part
);
1595 update_sd_second_half(gmx_int64_t step
,
1596 real
*dvdlambda
, /* the contribution to be added to the bonded interactions */
1597 const t_inputrec
*inputrec
, /* input record and box stuff */
1602 const t_commrec
*cr
,
1603 const gmx_multisim_t
*ms
,
1605 gmx_wallcycle_t wcycle
,
1607 gmx::Constraints
*constr
,
1615 if (inputrec
->eI
== eiSD1
)
1618 int homenr
= md
->homenr
;
1620 /* Cast delta_t from double to real to make the integrators faster.
1621 * The only reason for having delta_t double is to get accurate values
1622 * for t=delta_t*step when step is larger than float precision.
1623 * For integration dt the accuracy of real suffices, since with
1624 * integral += dt*integrand the increment is nearly always (much) smaller
1625 * than the integral (and the integrand has real precision).
1627 real dt
= inputrec
->delta_t
;
1629 wallcycle_start(wcycle
, ewcUPDATE
);
1631 nth
= gmx_omp_nthreads_get(emntUpdate
);
1633 #pragma omp parallel for num_threads(nth) schedule(static)
1634 for (th
= 0; th
< nth
; th
++)
1638 int start_th
, end_th
;
1639 getThreadAtomRange(nth
, th
, homenr
, &start_th
, &end_th
);
1641 doSDUpdateGeneral
<SDUpdate::FrictionAndNoiseOnly
>
1643 start_th
, end_th
, dt
,
1644 inputrec
->opts
.acc
, inputrec
->opts
.nFreeze
,
1645 md
->invmass
, md
->ptype
,
1646 md
->cFREEZE
, nullptr, md
->cTC
,
1647 as_rvec_array(state
->x
.data()), as_rvec_array(upd
->xp
.data()),
1648 as_rvec_array(state
->v
.data()), nullptr,
1649 step
, inputrec
->ld_seed
,
1650 DOMAINDECOMP(cr
) ? cr
->dd
->gatindex
: nullptr);
1652 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
1654 inc_nrnb(nrnb
, eNR_UPDATE
, homenr
);
1655 wallcycle_stop(wcycle
, ewcUPDATE
);
1658 /* Constrain the coordinates upd->xp for half a time step */
1659 wallcycle_start(wcycle
, ewcCONSTR
);
1660 constrain(nullptr, do_log
, do_ene
, constr
, idef
,
1661 inputrec
, cr
, ms
, step
, 1, 0.5, md
,
1662 as_rvec_array(state
->x
.data()), as_rvec_array(upd
->xp
.data()), nullptr,
1663 bMolPBC
, state
->box
,
1664 state
->lambda
[efptBONDED
], dvdlambda
,
1665 as_rvec_array(state
->v
.data()), nullptr, nrnb
, econqCoord
);
1667 wallcycle_stop(wcycle
, ewcCONSTR
);
1672 void finish_update(const t_inputrec
*inputrec
, /* input record and box stuff */
1677 gmx_wallcycle_t wcycle
,
1679 gmx::Constraints
*constr
)
1681 int homenr
= md
->homenr
;
1683 /* We must always unshift after updating coordinates; if we did not shake
1684 x was shifted in do_force */
1686 /* NOTE Currently we always integrate to a temporary buffer and
1687 * then copy the results back. */
1689 wallcycle_start_nocount(wcycle
, ewcUPDATE
);
1691 if (md
->cFREEZE
!= nullptr && constr
!= nullptr)
1693 /* If we have atoms that are frozen along some, but not all
1694 * dimensions, then any constraints will have moved them also along
1695 * the frozen dimensions. To freeze such degrees of freedom
1696 * we copy them back here to later copy them forward. It would
1697 * be more elegant and slightly more efficient to copies zero
1698 * times instead of twice, but the graph case below prevents this.
1700 const ivec
*nFreeze
= inputrec
->opts
.nFreeze
;
1701 bool partialFreezeAndConstraints
= false;
1702 for (int g
= 0; g
< inputrec
->opts
.ngfrz
; g
++)
1704 int numFreezeDim
= nFreeze
[g
][XX
] + nFreeze
[g
][YY
] + nFreeze
[g
][ZZ
];
1705 if (numFreezeDim
> 0 && numFreezeDim
< 3)
1707 partialFreezeAndConstraints
= true;
1710 if (partialFreezeAndConstraints
)
1712 for (int i
= 0; i
< homenr
; i
++)
1714 int g
= md
->cFREEZE
[i
];
1716 for (int d
= 0; d
< DIM
; d
++)
1720 upd
->xp
[i
][d
] = state
->x
[i
][d
];
1727 if (graph
&& (graph
->nnodes
> 0))
1729 unshift_x(graph
, state
->box
, as_rvec_array(state
->x
.data()), as_rvec_array(upd
->xp
.data()));
1730 if (TRICLINIC(state
->box
))
1732 inc_nrnb(nrnb
, eNR_SHIFTX
, 2*graph
->nnodes
);
1736 inc_nrnb(nrnb
, eNR_SHIFTX
, graph
->nnodes
);
1741 /* The copy is performance sensitive, so use a bare pointer */
1742 rvec
*xp
= as_rvec_array(upd
->xp
.data());
1743 #ifndef __clang_analyzer__
1744 // cppcheck-suppress unreadVariable
1745 int gmx_unused nth
= gmx_omp_nthreads_get(emntUpdate
);
1747 #pragma omp parallel for num_threads(nth) schedule(static)
1748 for (int i
= 0; i
< homenr
; i
++)
1750 // Trivial statement, does not throw
1751 copy_rvec(xp
[i
], state
->x
[i
]);
1754 wallcycle_stop(wcycle
, ewcUPDATE
);
1756 /* ############# END the update of velocities and positions ######### */
1759 void update_pcouple_after_coordinates(FILE *fplog
,
1761 const t_inputrec
*inputrec
,
1762 const t_mdatoms
*md
,
1763 const matrix pressure
,
1764 const matrix forceVirial
,
1765 const matrix constraintVirial
,
1766 const matrix parrinellorahmanMu
,
1772 int homenr
= md
->homenr
;
1774 /* Cast to real for faster code, no loss in precision (see comment above) */
1775 real dt
= inputrec
->delta_t
;
1778 /* now update boxes */
1779 switch (inputrec
->epc
)
1783 case (epcBERENDSEN
):
1784 if (isPressureCouplingStep(step
, inputrec
))
1786 real dtpc
= inputrec
->nstpcouple
*dt
;
1788 berendsen_pcoupl(fplog
, step
, inputrec
, dtpc
,
1789 pressure
, state
->box
,
1790 forceVirial
, constraintVirial
,
1791 mu
, &state
->baros_integral
);
1792 berendsen_pscale(inputrec
, mu
, state
->box
, state
->box_rel
,
1793 start
, homenr
, as_rvec_array(state
->x
.data()),
1797 case (epcPARRINELLORAHMAN
):
1798 if (isPressureCouplingStep(step
, inputrec
))
1800 /* The box velocities were updated in do_pr_pcoupl,
1801 * but we dont change the box vectors until we get here
1802 * since we need to be able to shift/unshift above.
1804 real dtpc
= inputrec
->nstpcouple
*dt
;
1805 for (int i
= 0; i
< DIM
; i
++)
1807 for (int m
= 0; m
<= i
; m
++)
1809 state
->box
[i
][m
] += dtpc
*state
->boxv
[i
][m
];
1812 preserve_box_shape(inputrec
, state
->box_rel
, state
->box
);
1814 /* Scale the coordinates */
1815 for (int n
= start
; n
< start
+ homenr
; n
++)
1817 tmvmul_ur0(parrinellorahmanMu
, state
->x
[n
], state
->x
[n
]);
1822 switch (inputrec
->epct
)
1824 case (epctISOTROPIC
):
1825 /* DIM * eta = ln V. so DIM*eta_new = DIM*eta_old + DIM*dt*veta =>
1826 ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) =>
1827 Side length scales as exp(veta*dt) */
1829 msmul(state
->box
, std::exp(state
->veta
*dt
), state
->box
);
1831 /* Relate veta to boxv. veta = d(eta)/dT = (1/DIM)*1/V dV/dT.
1832 o If we assume isotropic scaling, and box length scaling
1833 factor L, then V = L^DIM (det(M)). So dV/dt = DIM
1834 L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt. The
1835 determinant of B is L^DIM det(M), and the determinant
1836 of dB/dt is (dL/dT)^DIM det (M). veta will be
1837 (det(dB/dT)/det(B))^(1/3). Then since M =
1838 B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */
1840 msmul(state
->box
, state
->veta
, state
->boxv
);
1850 if (inputrecDeform(inputrec
))
1852 deform(upd
, start
, homenr
, as_rvec_array(state
->x
.data()), state
->box
, inputrec
, step
);
1856 void update_coords(gmx_int64_t step
,
1857 t_inputrec
*inputrec
, /* input record and box stuff */
1860 gmx::PaddedArrayRef
<gmx::RVec
> f
, /* forces on home particles */
1862 gmx_ekindata_t
*ekind
,
1866 const t_commrec
*cr
, /* these shouldn't be here -- need to think about it */
1867 gmx::Constraints
*constr
)
1869 gmx_bool bDoConstr
= (nullptr != constr
);
1871 /* Running the velocity half does nothing except for velocity verlet */
1872 if ((UpdatePart
== etrtVELOCITY1
|| UpdatePart
== etrtVELOCITY2
) &&
1873 !EI_VV(inputrec
->eI
))
1875 gmx_incons("update_coords called for velocity without VV integrator");
1878 int homenr
= md
->homenr
;
1880 /* Cast to real for faster code, no loss in precision (see comment above) */
1881 real dt
= inputrec
->delta_t
;
1883 /* We need to update the NMR restraint history when time averaging is used */
1884 if (state
->flags
& (1<<estDISRE_RM3TAV
))
1886 update_disres_history(fcd
, &state
->hist
);
1888 if (state
->flags
& (1<<estORIRE_DTAV
))
1890 update_orires_history(fcd
, &state
->hist
);
1893 /* ############# START The update of velocities and positions ######### */
1894 int nth
= gmx_omp_nthreads_get(emntUpdate
);
1896 #pragma omp parallel for num_threads(nth) schedule(static)
1897 for (int th
= 0; th
< nth
; th
++)
1901 int start_th
, end_th
;
1902 getThreadAtomRange(nth
, th
, homenr
, &start_th
, &end_th
);
1904 /* Strictly speaking, we would only need this check with SIMD
1905 * and for the actual SIMD width. But since the code currently
1906 * always adds padding for GMX_REAL_MAX_SIMD_WIDTH, we check that.
1908 size_t gmx_used_in_debug homenrSimdPadded
;
1909 homenrSimdPadded
= ((homenr
+ GMX_REAL_MAX_SIMD_WIDTH
- 1)/GMX_REAL_MAX_SIMD_WIDTH
)*GMX_REAL_MAX_SIMD_WIDTH
;
1910 GMX_ASSERT(state
->x
.size() >= homenrSimdPadded
, "state->x needs to be padded for SIMD access");
1911 GMX_ASSERT(upd
->xp
.size() >= homenrSimdPadded
, "upd->xp needs to be padded for SIMD access");
1912 GMX_ASSERT(state
->v
.size() >= homenrSimdPadded
, "state->v needs to be padded for SIMD access");
1913 GMX_ASSERT(f
.size() >= homenrSimdPadded
, "f needs to be padded for SIMD access");
1915 const rvec
*x_rvec
= as_rvec_array(state
->x
.data());
1916 rvec
*xp_rvec
= as_rvec_array(upd
->xp
.data());
1917 rvec
*v_rvec
= as_rvec_array(state
->v
.data());
1918 const rvec
*f_rvec
= as_rvec_array(f
.data());
1920 switch (inputrec
->eI
)
1923 do_update_md(start_th
, end_th
, step
, dt
,
1924 inputrec
, md
, ekind
, state
->box
,
1925 x_rvec
, xp_rvec
, v_rvec
, f_rvec
,
1926 state
->nosehoover_vxi
.data(), M
);
1931 // With constraints, the SD update is done in 2 parts
1932 doSDUpdateGeneral
<SDUpdate::ForcesOnly
>
1934 start_th
, end_th
, dt
,
1935 inputrec
->opts
.acc
, inputrec
->opts
.nFreeze
,
1936 md
->invmass
, md
->ptype
,
1937 md
->cFREEZE
, md
->cACC
, nullptr,
1938 x_rvec
, xp_rvec
, v_rvec
, f_rvec
,
1939 step
, inputrec
->ld_seed
, nullptr);
1943 doSDUpdateGeneral
<SDUpdate::Combined
>
1945 start_th
, end_th
, dt
,
1946 inputrec
->opts
.acc
, inputrec
->opts
.nFreeze
,
1947 md
->invmass
, md
->ptype
,
1948 md
->cFREEZE
, md
->cACC
, md
->cTC
,
1949 x_rvec
, xp_rvec
, v_rvec
, f_rvec
,
1950 step
, inputrec
->ld_seed
,
1951 DOMAINDECOMP(cr
) ? cr
->dd
->gatindex
: nullptr);
1955 do_update_bd(start_th
, end_th
, dt
,
1956 inputrec
->opts
.nFreeze
, md
->invmass
, md
->ptype
,
1957 md
->cFREEZE
, md
->cTC
,
1958 x_rvec
, xp_rvec
, v_rvec
, f_rvec
,
1961 step
, inputrec
->ld_seed
, DOMAINDECOMP(cr
) ? cr
->dd
->gatindex
: nullptr);
1966 gmx_bool bExtended
= (inputrec
->etc
== etcNOSEHOOVER
||
1967 inputrec
->epc
== epcPARRINELLORAHMAN
||
1968 inputrec
->epc
== epcMTTK
);
1970 /* assuming barostat coupled to group 0 */
1971 real alpha
= 1.0 + DIM
/static_cast<real
>(inputrec
->opts
.nrdf
[0]);
1976 do_update_vv_vel(start_th
, end_th
, dt
,
1977 inputrec
->opts
.acc
, inputrec
->opts
.nFreeze
,
1978 md
->invmass
, md
->ptype
,
1979 md
->cFREEZE
, md
->cACC
,
1981 bExtended
, state
->veta
, alpha
);
1984 do_update_vv_pos(start_th
, end_th
, dt
,
1985 inputrec
->opts
.nFreeze
,
1986 md
->ptype
, md
->cFREEZE
,
1987 x_rvec
, xp_rvec
, v_rvec
,
1988 bExtended
, state
->veta
);
1994 gmx_fatal(FARGS
, "Don't know how to update coordinates");
1998 GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR
;
2003 extern gmx_bool
update_randomize_velocities(t_inputrec
*ir
, gmx_int64_t step
, const t_commrec
*cr
,
2004 t_mdatoms
*md
, t_state
*state
, gmx_update_t
*upd
, gmx::Constraints
*constr
)
2007 real rate
= (ir
->delta_t
)/ir
->opts
.tau_t
[0];
2009 if (ir
->etc
== etcANDERSEN
&& constr
!= nullptr)
2011 /* Currently, Andersen thermostat does not support constrained
2012 systems. Functionality exists in the andersen_tcoupl
2013 function in GROMACS 4.5.7 to allow this combination. That
2014 code could be ported to the current random-number
2015 generation approach, but has not yet been done because of
2016 lack of time and resources. */
2017 gmx_fatal(FARGS
, "Normal Andersen is currently not supported with constraints, use massive Andersen instead");
2020 /* proceed with andersen if 1) it's fixed probability per
2021 particle andersen or 2) it's massive andersen and it's tau_t/dt */
2022 if ((ir
->etc
== etcANDERSEN
) || do_per_step(step
, static_cast<int>(1.0/rate
+ 0.5)))
2024 andersen_tcoupl(ir
, step
, cr
, md
, state
, rate
,
2025 upd
->sd
->randomize_group
, upd
->sd
->boltzfac
);