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,2019,2020, 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.
44 #include "gromacs/domdec/domdec_struct.h"
45 #include "gromacs/gmxlib/nrnb.h"
46 #include "gromacs/math/functions.h"
47 #include "gromacs/math/invertmatrix.h"
48 #include "gromacs/math/units.h"
49 #include "gromacs/math/vec.h"
50 #include "gromacs/math/vecdump.h"
51 #include "gromacs/mdlib/expanded.h"
52 #include "gromacs/mdlib/gmx_omp_nthreads.h"
53 #include "gromacs/mdlib/stat.h"
54 #include "gromacs/mdlib/update.h"
55 #include "gromacs/mdtypes/commrec.h"
56 #include "gromacs/mdtypes/enerdata.h"
57 #include "gromacs/mdtypes/group.h"
58 #include "gromacs/mdtypes/inputrec.h"
59 #include "gromacs/mdtypes/md_enums.h"
60 #include "gromacs/mdtypes/mdatom.h"
61 #include "gromacs/mdtypes/state.h"
62 #include "gromacs/pbcutil/boxutilities.h"
63 #include "gromacs/pbcutil/pbc.h"
64 #include "gromacs/random/gammadistribution.h"
65 #include "gromacs/random/normaldistribution.h"
66 #include "gromacs/random/tabulatednormaldistribution.h"
67 #include "gromacs/random/threefry.h"
68 #include "gromacs/random/uniformrealdistribution.h"
69 #include "gromacs/utility/cstringutil.h"
70 #include "gromacs/utility/fatalerror.h"
71 #include "gromacs/utility/pleasecite.h"
72 #include "gromacs/utility/smalloc.h"
74 #define NTROTTERPARTS 3
76 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
78 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
79 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
81 #define MAX_SUZUKI_YOSHIDA_NUM 5
82 #define SUZUKI_YOSHIDA_NUM 5
84 static const double sy_const_1[] = { 1. };
85 static const double sy_const_3[] = { 0.828981543588751, -0.657963087177502, 0.828981543588751 };
86 static const double sy_const_5[] = { 0.2967324292201065, 0.2967324292201065, -0.186929716880426,
87 0.2967324292201065, 0.2967324292201065 };
89 static const double* sy_const[] = { nullptr, sy_const_1, nullptr, sy_const_3, nullptr, sy_const_5 };
92 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
96 {0.828981543588751,-0.657963087177502,0.828981543588751},
98 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
101 /* these integration routines are only referenced inside this file */
102 static void NHC_trotter(const t_grpopts* opts,
104 const gmx_ekindata_t* ekind,
110 const t_extmass* MassQ,
111 gmx_bool bEkinAveVel)
114 /* general routine for both barostat and thermostat nose hoover chains */
117 double Ekin, Efac, reft, kT, nd;
122 int mstepsi, mstepsj;
123 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
124 int nh = opts->nhchainlength;
127 mstepsi = mstepsj = ns;
129 /* if scalefac is NULL, we are doing the NHC of the barostat */
132 if (scalefac == nullptr)
137 for (i = 0; i < nvar; i++)
140 /* make it easier to iterate by selecting
141 out the sub-array that corresponds to this T group */
145 gmx::ArrayRef<const double> iQinv;
148 iQinv = gmx::arrayRefFromArray(&MassQ->QPinv[i * nh], nh);
149 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
150 reft = std::max<real>(0, opts->ref_t[0]);
151 Ekin = gmx::square(*veta) / MassQ->Winv;
155 iQinv = gmx::arrayRefFromArray(&MassQ->Qinv[i * nh], nh);
156 const t_grp_tcstat* tcstat = &ekind->tcstat[i];
158 reft = std::max<real>(0, opts->ref_t[i]);
161 Ekin = 2 * trace(tcstat->ekinf) * tcstat->ekinscalef_nhc;
165 Ekin = 2 * trace(tcstat->ekinh) * tcstat->ekinscaleh_nhc;
170 for (mi = 0; mi < mstepsi; mi++)
172 for (mj = 0; mj < mstepsj; mj++)
174 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
175 dt = sy_const[ns][mj] * dtfull / mstepsi;
177 /* compute the thermal forces */
178 GQ[0] = iQinv[0] * (Ekin - nd * kT);
180 for (j = 0; j < nh - 1; j++)
182 if (iQinv[j + 1] > 0)
184 /* we actually don't need to update here if we save the
185 state of the GQ, but it's easier to just recompute*/
186 GQ[j + 1] = iQinv[j + 1] * ((gmx::square(ivxi[j]) / iQinv[j]) - kT);
194 ivxi[nh - 1] += 0.25 * dt * GQ[nh - 1];
195 for (j = nh - 1; j > 0; j--)
197 Efac = exp(-0.125 * dt * ivxi[j]);
198 ivxi[j - 1] = Efac * (ivxi[j - 1] * Efac + 0.25 * dt * GQ[j - 1]);
201 Efac = exp(-0.5 * dt * ivxi[0]);
210 Ekin *= (Efac * Efac);
212 /* Issue - if the KE is an average of the last and the current temperatures, then we
213 might not be able to scale the kinetic energy directly with this factor. Might
214 take more bookkeeping -- have to think about this a bit more . . . */
216 GQ[0] = iQinv[0] * (Ekin - nd * kT);
218 /* update thermostat positions */
219 for (j = 0; j < nh; j++)
221 ixi[j] += 0.5 * dt * ivxi[j];
224 for (j = 0; j < nh - 1; j++)
226 Efac = exp(-0.125 * dt * ivxi[j + 1]);
227 ivxi[j] = Efac * (ivxi[j] * Efac + 0.25 * dt * GQ[j]);
228 if (iQinv[j + 1] > 0)
230 GQ[j + 1] = iQinv[j + 1] * ((gmx::square(ivxi[j]) / iQinv[j]) - kT);
237 ivxi[nh - 1] += 0.25 * dt * GQ[nh - 1];
244 static void boxv_trotter(const t_inputrec* ir,
248 const gmx_ekindata_t* ekind,
251 const t_extmass* MassQ)
258 tensor ekinmod, localpres;
260 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
261 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
264 if (ir->epct == epctSEMIISOTROPIC)
273 /* eta is in pure units. veta is in units of ps^-1. GW is in
274 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
275 taken to use only RATIOS of eta in updating the volume. */
277 /* we take the partial pressure tensors, modify the
278 kinetic energy tensor, and recovert to pressure */
280 if (ir->opts.nrdf[0] == 0)
282 gmx_fatal(FARGS, "Barostat is coupled to a T-group with no degrees of freedom\n");
284 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
285 alpha = 1.0 + DIM / (static_cast<double>(ir->opts.nrdf[0]));
286 alpha *= ekind->tcstat[0].ekinscalef_nhc;
287 msmul(ekind->ekin, alpha, ekinmod);
288 /* for now, we use Elr = 0, because if you want to get it right, you
289 really should be using PME. Maybe print a warning? */
291 pscal = calc_pres(ir->ePBC, nwall, box, ekinmod, vir, localpres) + pcorr;
294 GW = (vol * (MassQ->Winv / PRESFAC)) * (DIM * pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
296 *veta += 0.5 * dt * GW;
300 * This file implements temperature and pressure coupling algorithms:
301 * For now only the Weak coupling and the modified weak coupling.
303 * Furthermore computation of pressure and temperature is done here
307 real calc_pres(int ePBC, int nwall, const matrix box, const tensor ekin, const tensor vir, tensor pres)
312 if (ePBC == epbcNONE || (ePBC == epbcXY && nwall != 2))
318 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
319 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
323 fac = PRESFAC * 2.0 / det(box);
324 for (n = 0; (n < DIM); n++)
326 for (m = 0; (m < DIM); m++)
328 pres[n][m] = (ekin[n][m] - vir[n][m]) * fac;
334 pr_rvecs(debug, 0, "PC: pres", pres, DIM);
335 pr_rvecs(debug, 0, "PC: ekin", ekin, DIM);
336 pr_rvecs(debug, 0, "PC: vir ", vir, DIM);
337 pr_rvecs(debug, 0, "PC: box ", box, DIM);
340 return trace(pres) / DIM;
343 real calc_temp(real ekin, real nrdf)
347 return (2.0 * ekin) / (nrdf * BOLTZ);
355 /*! \brief Sets 1/mass for Parrinello-Rahman in wInv; NOTE: PRESFAC is not included, so not in GROMACS units! */
356 static void calcParrinelloRahmanInvMass(const t_inputrec* ir, const matrix box, tensor wInv)
360 /* TODO: See if we can make the mass independent of the box size */
361 maxBoxLength = std::max(box[XX][XX], box[YY][YY]);
362 maxBoxLength = std::max(maxBoxLength, box[ZZ][ZZ]);
364 for (int d = 0; d < DIM; d++)
366 for (int n = 0; n < DIM; n++)
368 wInv[d][n] = (4 * M_PI * M_PI * ir->compress[d][n]) / (3 * ir->tau_p * ir->tau_p * maxBoxLength);
373 void parrinellorahman_pcoupl(FILE* fplog,
375 const t_inputrec* ir,
385 /* This doesn't do any coordinate updating. It just
386 * integrates the box vector equations from the calculated
387 * acceleration due to pressure difference. We also compute
388 * the tensor M which is used in update to couple the particle
389 * coordinates to the box vectors.
391 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
394 * M_nk = (h') * (h' * h + h' h) * h
396 * with the dots denoting time derivatives and h is the transformation from
397 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
398 * This also goes for the pressure and M tensors - they are transposed relative
399 * to ours. Our equation thus becomes:
402 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
404 * where b is the gromacs box matrix.
405 * Our box accelerations are given by
407 * b = vol/W inv(box') * (P-ref_P) (=h')
410 real vol = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
411 real atot, arel, change;
412 tensor invbox, pdiff, t1, t2;
414 gmx::invertBoxMatrix(box, invbox);
418 /* Note that PRESFAC does not occur here.
419 * The pressure and compressibility always occur as a product,
420 * therefore the pressure unit drops out.
423 calcParrinelloRahmanInvMass(ir, box, winv);
425 m_sub(pres, ir->ref_p, pdiff);
427 if (ir->epct == epctSURFACETENSION)
429 /* Unlike Berendsen coupling it might not be trivial to include a z
430 * pressure correction here? On the other hand we don't scale the
431 * box momentarily, but change accelerations, so it might not be crucial.
433 real xy_pressure = 0.5 * (pres[XX][XX] + pres[YY][YY]);
434 for (int d = 0; d < ZZ; d++)
436 pdiff[d][d] = (xy_pressure - (pres[ZZ][ZZ] - ir->ref_p[d][d] / box[d][d]));
440 tmmul(invbox, pdiff, t1);
441 /* Move the off-diagonal elements of the 'force' to one side to ensure
442 * that we obey the box constraints.
444 for (int d = 0; d < DIM; d++)
446 for (int n = 0; n < d; n++)
448 t1[d][n] += t1[n][d];
455 case epctANISOTROPIC:
456 for (int d = 0; d < DIM; d++)
458 for (int n = 0; n <= d; n++)
460 t1[d][n] *= winv[d][n] * vol;
465 /* calculate total volume acceleration */
466 atot = box[XX][XX] * box[YY][YY] * t1[ZZ][ZZ] + box[XX][XX] * t1[YY][YY] * box[ZZ][ZZ]
467 + t1[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
468 arel = atot / (3 * vol);
469 /* set all RELATIVE box accelerations equal, and maintain total V
471 for (int d = 0; d < DIM; d++)
473 for (int n = 0; n <= d; n++)
475 t1[d][n] = winv[0][0] * vol * arel * box[d][n];
479 case epctSEMIISOTROPIC:
480 case epctSURFACETENSION:
481 /* Note the correction to pdiff above for surftens. coupling */
483 /* calculate total XY volume acceleration */
484 atot = box[XX][XX] * t1[YY][YY] + t1[XX][XX] * box[YY][YY];
485 arel = atot / (2 * box[XX][XX] * box[YY][YY]);
486 /* set RELATIVE XY box accelerations equal, and maintain total V
487 * change speed. Dont change the third box vector accelerations */
488 for (int d = 0; d < ZZ; d++)
490 for (int n = 0; n <= d; n++)
492 t1[d][n] = winv[d][n] * vol * arel * box[d][n];
495 for (int n = 0; n < DIM; n++)
497 t1[ZZ][n] *= winv[ZZ][n] * vol;
502 "Parrinello-Rahman pressure coupling type %s "
503 "not supported yet\n",
504 EPCOUPLTYPETYPE(ir->epct));
508 for (int d = 0; d < DIM; d++)
510 for (int n = 0; n <= d; n++)
512 boxv[d][n] += dt * t1[d][n];
514 /* We do NOT update the box vectors themselves here, since
515 * we need them for shifting later. It is instead done last
516 * in the update() routine.
519 /* Calculate the change relative to diagonal elements-
520 since it's perfectly ok for the off-diagonal ones to
521 be zero it doesn't make sense to check the change relative
525 change = std::fabs(dt * boxv[d][n] / box[d][d]);
527 if (change > maxchange)
534 if (maxchange > 0.01 && fplog)
538 "\nStep %s Warning: Pressure scaling more than 1%%. "
539 "This may mean your system\n is not yet equilibrated. "
540 "Use of Parrinello-Rahman pressure coupling during\n"
541 "equilibration can lead to simulation instability, "
542 "and is discouraged.\n",
543 gmx_step_str(step, buf));
547 preserve_box_shape(ir, box_rel, boxv);
549 mtmul(boxv, box, t1); /* t1=boxv * b' */
550 mmul(invbox, t1, t2);
551 mtmul(t2, invbox, M);
553 /* Determine the scaling matrix mu for the coordinates */
554 for (int d = 0; d < DIM; d++)
556 for (int n = 0; n <= d; n++)
558 t1[d][n] = box[d][n] + dt * boxv[d][n];
561 preserve_box_shape(ir, box_rel, t1);
562 /* t1 is the box at t+dt, determine mu as the relative change */
563 mmul_ur0(invbox, t1, mu);
566 void berendsen_pcoupl(FILE* fplog,
568 const t_inputrec* ir,
572 const matrix force_vir,
573 const matrix constraint_vir,
575 double* baros_integral)
578 real scalar_pressure, xy_pressure, p_corr_z;
582 * Calculate the scaling matrix mu
586 for (d = 0; d < DIM; d++)
588 scalar_pressure += pres[d][d] / DIM;
591 xy_pressure += pres[d][d] / (DIM - 1);
594 /* Pressure is now in bar, everywhere. */
595 #define factor(d, m) (ir->compress[d][m] * dt / ir->tau_p)
597 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
598 * necessary for triclinic scaling
604 for (d = 0; d < DIM; d++)
606 mu[d][d] = 1.0 - factor(d, d) * (ir->ref_p[d][d] - scalar_pressure) / DIM;
609 case epctSEMIISOTROPIC:
610 for (d = 0; d < ZZ; d++)
612 mu[d][d] = 1.0 - factor(d, d) * (ir->ref_p[d][d] - xy_pressure) / DIM;
614 mu[ZZ][ZZ] = 1.0 - factor(ZZ, ZZ) * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM;
616 case epctANISOTROPIC:
617 for (d = 0; d < DIM; d++)
619 for (n = 0; n < DIM; n++)
621 mu[d][n] = (d == n ? 1.0 : 0.0) - factor(d, n) * (ir->ref_p[d][n] - pres[d][n]) / DIM;
625 case epctSURFACETENSION:
626 /* ir->ref_p[0/1] is the reference surface-tension times *
627 * the number of surfaces */
628 if (ir->compress[ZZ][ZZ] != 0.0F)
630 p_corr_z = dt / ir->tau_p * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
634 /* when the compressibity is zero, set the pressure correction *
635 * in the z-direction to zero to get the correct surface tension */
638 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ] * p_corr_z;
639 for (d = 0; d < DIM - 1; d++)
643 * (ir->ref_p[d][d] / (mu[ZZ][ZZ] * box[ZZ][ZZ])
644 - (pres[ZZ][ZZ] + p_corr_z - xy_pressure))
649 gmx_fatal(FARGS, "Berendsen pressure coupling type %s not supported yet\n",
650 EPCOUPLTYPETYPE(ir->epct));
652 /* To fullfill the orientation restrictions on triclinic boxes
653 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
654 * the other elements of mu to first order.
656 mu[YY][XX] += mu[XX][YY];
657 mu[ZZ][XX] += mu[XX][ZZ];
658 mu[ZZ][YY] += mu[YY][ZZ];
663 /* Keep track of the work the barostat applies on the system.
664 * Without constraints force_vir tells us how Epot changes when scaling.
665 * With constraints constraint_vir gives us the constraint contribution
666 * to both Epot and Ekin. Although we are not scaling velocities, scaling
667 * the coordinates leads to scaling of distances involved in constraints.
668 * This in turn changes the angular momentum (even if the constrained
669 * distances are corrected at the next step). The kinetic component
670 * of the constraint virial captures the angular momentum change.
672 for (int d = 0; d < DIM; d++)
674 for (int n = 0; n <= d; n++)
677 2 * (mu[d][n] - (n == d ? 1 : 0)) * (force_vir[d][n] + constraint_vir[d][n]);
683 pr_rvecs(debug, 0, "PC: pres ", pres, 3);
684 pr_rvecs(debug, 0, "PC: mu ", mu, 3);
687 if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 || mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01
688 || mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
692 "\nStep %s Warning: pressure scaling more than 1%%, "
694 gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
697 fprintf(fplog, "%s", buf);
699 fprintf(stderr, "%s", buf);
703 void berendsen_pscale(const t_inputrec* ir,
710 const unsigned short cFREEZE[],
712 const bool scaleCoordinates)
714 ivec* nFreeze = ir->opts.nFreeze;
716 int nthreads gmx_unused;
718 #ifndef __clang_analyzer__
719 nthreads = gmx_omp_nthreads_get(emntUpdate);
722 /* Scale the positions */
723 if (scaleCoordinates)
725 #pragma omp parallel for num_threads(nthreads) schedule(static)
726 for (int n = start; n < start + nr_atoms; n++)
728 // Trivial OpenMP region that does not throw
731 if (cFREEZE == nullptr)
742 x[n][XX] = mu[XX][XX] * x[n][XX] + mu[YY][XX] * x[n][YY] + mu[ZZ][XX] * x[n][ZZ];
746 x[n][YY] = mu[YY][YY] * x[n][YY] + mu[ZZ][YY] * x[n][ZZ];
750 x[n][ZZ] = mu[ZZ][ZZ] * x[n][ZZ];
754 /* compute final boxlengths */
755 for (d = 0; d < DIM; d++)
757 box[d][XX] = mu[XX][XX] * box[d][XX] + mu[YY][XX] * box[d][YY] + mu[ZZ][XX] * box[d][ZZ];
758 box[d][YY] = mu[YY][YY] * box[d][YY] + mu[ZZ][YY] * box[d][ZZ];
759 box[d][ZZ] = mu[ZZ][ZZ] * box[d][ZZ];
762 preserve_box_shape(ir, box_rel, box);
764 /* (un)shifting should NOT be done after this,
765 * since the box vectors might have changed
767 inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
770 void berendsen_tcoupl(const t_inputrec* ir, gmx_ekindata_t* ekind, real dt, std::vector<double>& therm_integral)
772 const t_grpopts* opts = &ir->opts;
774 for (int i = 0; (i < opts->ngtc); i++)
780 Ek = trace(ekind->tcstat[i].ekinf);
781 T = ekind->tcstat[i].T;
785 Ek = trace(ekind->tcstat[i].ekinh);
786 T = ekind->tcstat[i].Th;
789 if ((opts->tau_t[i] > 0) && (T > 0.0))
791 real reft = std::max<real>(0, opts->ref_t[i]);
792 real lll = std::sqrt(1.0 + (dt / opts->tau_t[i]) * (reft / T - 1.0));
793 ekind->tcstat[i].lambda = std::max<real>(std::min<real>(lll, 1.25), 0.8);
797 ekind->tcstat[i].lambda = 1.0;
800 /* Keep track of the amount of energy we are adding to the system */
801 therm_integral[i] -= (gmx::square(ekind->tcstat[i].lambda) - 1) * Ek;
805 fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n", i, T, ekind->tcstat[i].lambda);
810 void andersen_tcoupl(const t_inputrec* ir,
814 gmx::ArrayRef<gmx::RVec> v,
816 const std::vector<bool>& randomize,
817 gmx::ArrayRef<const real> boltzfac)
819 const int* gatindex = (DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
822 gmx::ThreeFry2x64<0> rng(ir->andersen_seed, gmx::RandomDomain::Thermostat);
823 gmx::UniformRealDistribution<real> uniformDist;
824 gmx::TabulatedNormalDistribution<real, 14> normalDist;
826 /* randomize the velocities of the selected particles */
828 for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */
830 int ng = gatindex ? gatindex[i] : i;
833 rng.restart(step, ng);
837 gc = md->cTC[i]; /* assign the atom to a temperature group if there are more than one */
841 if (ir->etc == etcANDERSENMASSIVE)
843 /* Randomize particle always */
848 /* Randomize particle probabilistically */
850 bRandomize = uniformDist(rng) < rate;
857 scal = std::sqrt(boltzfac[gc] * md->invmass[i]);
861 for (d = 0; d < DIM; d++)
863 v[i][d] = scal * normalDist(rng);
871 void nosehoover_tcoupl(const t_grpopts* opts,
872 const gmx_ekindata_t* ekind,
876 const t_extmass* MassQ)
881 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
883 for (i = 0; (i < opts->ngtc); i++)
885 reft = std::max<real>(0, opts->ref_t[i]);
887 vxi[i] += dt * MassQ->Qinv[i] * (ekind->tcstat[i].Th - reft);
888 xi[i] += dt * (oldvxi + vxi[i]) * 0.5;
892 void trotter_update(const t_inputrec* ir,
894 gmx_ekindata_t* ekind,
895 const gmx_enerdata_t* enerd,
899 const t_extmass* MassQ,
900 gmx::ArrayRef<std::vector<int>> trotter_seqlist,
904 int n, i, d, ngtc, gc = 0, t;
905 t_grp_tcstat* tcstat;
906 const t_grpopts* opts;
909 double * scalefac, dtc;
910 rvec sumv = { 0, 0, 0 };
913 if (trotter_seqno <= ettTSEQ2)
915 step_eff = step - 1; /* the velocity verlet calls are actually out of order -- the first
916 half step is actually the last half step from the previous step.
917 Thus the first half step actually corresponds to the n-1 step*/
924 bCouple = (ir->nsttcouple == 1 || do_per_step(step_eff + ir->nsttcouple, ir->nsttcouple));
926 const gmx::ArrayRef<const int> trotter_seq = trotter_seqlist[trotter_seqno];
928 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
932 dtc = ir->nsttcouple * ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
933 opts = &(ir->opts); /* just for ease of referencing */
936 snew(scalefac, opts->ngtc);
937 for (i = 0; i < ngtc; i++)
941 /* execute the series of trotter updates specified in the trotterpart array */
943 for (i = 0; i < NTROTTERPARTS; i++)
945 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
946 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2)
947 || (trotter_seq[i] == etrtNHC2))
956 auto v = makeArrayRef(state->v);
957 switch (trotter_seq[i])
961 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir,
962 enerd->term[F_PDISPCORR], MassQ);
966 NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi.data(),
967 state->nhpres_vxi.data(), nullptr, &(state->veta), MassQ, FALSE);
971 NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi.data(),
972 state->nosehoover_vxi.data(), scalefac, nullptr, MassQ, (ir->eI == eiVV));
973 /* need to rescale the kinetic energies and velocities here. Could
974 scale the velocities later, but we need them scaled in order to
975 produce the correct outputs, so we'll scale them here. */
977 for (t = 0; t < ngtc; t++)
979 tcstat = &ekind->tcstat[t];
980 tcstat->vscale_nhc = scalefac[t];
981 tcstat->ekinscaleh_nhc *= (scalefac[t] * scalefac[t]);
982 tcstat->ekinscalef_nhc *= (scalefac[t] * scalefac[t]);
984 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
985 /* but do we actually need the total? */
987 /* modify the velocities as well */
988 for (n = 0; n < md->homenr; n++)
990 if (md->cTC) /* does this conditional need to be here? is this always true?*/
994 for (d = 0; d < DIM; d++)
996 v[n][d] *= scalefac[gc];
1001 for (d = 0; d < DIM; d++)
1003 sumv[d] += (v[n][d]) / md->invmass[n];
1011 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
1016 extern void init_npt_masses(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool bInit)
1018 int n, i, j, d, ngtc, nh;
1019 const t_grpopts* opts;
1020 real reft, kT, ndj, nd;
1022 opts = &(ir->opts); /* just for ease of referencing */
1023 ngtc = ir->opts.ngtc;
1024 nh = state->nhchainlength;
1030 MassQ->Qinv.resize(ngtc);
1032 for (i = 0; (i < ngtc); i++)
1034 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1036 MassQ->Qinv[i] = 1.0 / (gmx::square(opts->tau_t[i] / M_2PI) * opts->ref_t[i]);
1040 MassQ->Qinv[i] = 0.0;
1044 else if (EI_VV(ir->eI))
1046 /* Set pressure variables */
1050 if (state->vol0 == 0)
1052 state->vol0 = det(state->box);
1053 /* because we start by defining a fixed
1054 compressibility, we need the volume at this
1055 compressibility to solve the problem. */
1059 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1060 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1061 MassQ->Winv = (PRESFAC * trace(ir->compress) * BOLTZ * opts->ref_t[0])
1062 / (DIM * state->vol0 * gmx::square(ir->tau_p / M_2PI));
1063 /* An alternate mass definition, from Tuckerman et al. */
1064 /* MassQ->Winv = 1.0/(gmx::square(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1065 for (d = 0; d < DIM; d++)
1067 for (n = 0; n < DIM; n++)
1069 MassQ->Winvm[d][n] =
1070 PRESFAC * ir->compress[d][n] / (state->vol0 * gmx::square(ir->tau_p / M_2PI));
1071 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1072 before using MTTK for anisotropic states.*/
1075 /* Allocate space for thermostat variables */
1078 MassQ->Qinv.resize(ngtc * nh);
1081 /* now, set temperature variables */
1082 for (i = 0; i < ngtc; i++)
1084 if (opts->tau_t[i] > 0 && opts->ref_t[i] > 0 && opts->nrdf[i] > 0)
1086 reft = std::max<real>(0, opts->ref_t[i]);
1089 for (j = 0; j < nh; j++)
1099 MassQ->Qinv[i * nh + j] = 1.0 / (gmx::square(opts->tau_t[i] / M_2PI) * ndj * kT);
1104 for (j = 0; j < nh; j++)
1106 MassQ->Qinv[i * nh + j] = 0.0;
1113 std::array<std::vector<int>, ettTSEQMAX>
1114 init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool bTrotter)
1116 int i, j, nnhpres, nh;
1117 const t_grpopts* opts;
1118 real bmass, qmass, reft, kT;
1120 opts = &(ir->opts); /* just for ease of referencing */
1121 nnhpres = state->nnhpres;
1122 nh = state->nhchainlength;
1124 if (EI_VV(ir->eI) && (ir->epc == epcMTTK) && (ir->etc != etcNOSEHOOVER))
1127 "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
1130 init_npt_masses(ir, state, MassQ, TRUE);
1132 /* first, initialize clear all the trotter calls */
1133 std::array<std::vector<int>, ettTSEQMAX> trotter_seq;
1134 for (i = 0; i < ettTSEQMAX; i++)
1136 trotter_seq[i].resize(NTROTTERPARTS, etrtNONE);
1137 trotter_seq[i][0] = etrtSKIPALL;
1142 /* no trotter calls, so we never use the values in the array.
1143 * We access them (so we need to define them, but ignore
1149 /* compute the kinetic energy by using the half step velocities or
1150 * the kinetic energies, depending on the order of the trotter calls */
1154 if (inputrecNptTrotter(ir))
1156 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1157 We start with the initial one. */
1158 /* first, a round that estimates veta. */
1159 trotter_seq[0][0] = etrtBAROV;
1161 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1163 /* The first half trotter update */
1164 trotter_seq[2][0] = etrtBAROV;
1165 trotter_seq[2][1] = etrtNHC;
1166 trotter_seq[2][2] = etrtBARONHC;
1168 /* The second half trotter update */
1169 trotter_seq[3][0] = etrtBARONHC;
1170 trotter_seq[3][1] = etrtNHC;
1171 trotter_seq[3][2] = etrtBAROV;
1173 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1175 else if (inputrecNvtTrotter(ir))
1177 /* This is the easy version - there are only two calls, both the same.
1178 Otherwise, even easier -- no calls */
1179 trotter_seq[2][0] = etrtNHC;
1180 trotter_seq[3][0] = etrtNHC;
1182 else if (inputrecNphTrotter(ir))
1184 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1185 We start with the initial one. */
1186 /* first, a round that estimates veta. */
1187 trotter_seq[0][0] = etrtBAROV;
1189 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1191 /* The first half trotter update */
1192 trotter_seq[2][0] = etrtBAROV;
1193 trotter_seq[2][1] = etrtBARONHC;
1195 /* The second half trotter update */
1196 trotter_seq[3][0] = etrtBARONHC;
1197 trotter_seq[3][1] = etrtBAROV;
1199 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1202 else if (ir->eI == eiVVAK)
1204 if (inputrecNptTrotter(ir))
1206 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1207 We start with the initial one. */
1208 /* first, a round that estimates veta. */
1209 trotter_seq[0][0] = etrtBAROV;
1211 /* The first half trotter update, part 1 -- double update, because it commutes */
1212 trotter_seq[1][0] = etrtNHC;
1214 /* The first half trotter update, part 2 */
1215 trotter_seq[2][0] = etrtBAROV;
1216 trotter_seq[2][1] = etrtBARONHC;
1218 /* The second half trotter update, part 1 */
1219 trotter_seq[3][0] = etrtBARONHC;
1220 trotter_seq[3][1] = etrtBAROV;
1222 /* The second half trotter update */
1223 trotter_seq[4][0] = etrtNHC;
1225 else if (inputrecNvtTrotter(ir))
1227 /* This is the easy version - there is only one call, both the same.
1228 Otherwise, even easier -- no calls */
1229 trotter_seq[1][0] = etrtNHC;
1230 trotter_seq[4][0] = etrtNHC;
1232 else if (inputrecNphTrotter(ir))
1234 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1235 We start with the initial one. */
1236 /* first, a round that estimates veta. */
1237 trotter_seq[0][0] = etrtBAROV;
1239 /* The first half trotter update, part 1 -- leave zero */
1240 trotter_seq[1][0] = etrtNHC;
1242 /* The first half trotter update, part 2 */
1243 trotter_seq[2][0] = etrtBAROV;
1244 trotter_seq[2][1] = etrtBARONHC;
1246 /* The second half trotter update, part 1 */
1247 trotter_seq[3][0] = etrtBARONHC;
1248 trotter_seq[3][1] = etrtBAROV;
1250 /* The second half trotter update -- blank for now */
1257 default: bmass = DIM * DIM; /* recommended mass parameters for isotropic barostat */
1260 MassQ->QPinv.resize(nnhpres * opts->nhchainlength);
1262 /* barostat temperature */
1263 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1265 reft = std::max<real>(0, opts->ref_t[0]);
1267 for (i = 0; i < nnhpres; i++)
1269 for (j = 0; j < nh; j++)
1279 MassQ->QPinv[i * opts->nhchainlength + j] =
1280 1.0 / (gmx::square(opts->tau_t[0] / M_2PI) * qmass * kT);
1286 for (i = 0; i < nnhpres; i++)
1288 for (j = 0; j < nh; j++)
1290 MassQ->QPinv[i * nh + j] = 0.0;
1297 static real energyNoseHoover(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1301 int nh = state->nhchainlength;
1303 for (int i = 0; i < ir->opts.ngtc; i++)
1305 const double* ixi = &state->nosehoover_xi[i * nh];
1306 const double* ivxi = &state->nosehoover_vxi[i * nh];
1307 const double* iQinv = &(MassQ->Qinv[i * nh]);
1309 real nd = ir->opts.nrdf[i];
1310 real reft = std::max<real>(ir->opts.ref_t[i], 0);
1311 real kT = BOLTZ * reft;
1315 if (inputrecNvtTrotter(ir) || inputrecNptTrotter(ir))
1317 /* contribution from the thermal momenta of the NH chain */
1318 for (int j = 0; j < nh; j++)
1322 energy += 0.5 * gmx::square(ivxi[j]) / iQinv[j];
1323 /* contribution from the thermal variable of the NH chain */
1333 energy += ndj * ixi[j] * kT;
1337 else /* Other non Trotter temperature NH control -- no chains yet. */
1339 energy += 0.5 * BOLTZ * nd * gmx::square(ivxi[0]) / iQinv[0];
1340 energy += nd * ixi[0] * kT;
1348 /* Returns the energy from the barostat thermostat chain */
1349 static real energyPressureMTTK(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1353 int nh = state->nhchainlength;
1355 for (int i = 0; i < state->nnhpres; i++)
1357 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1358 real reft = std::max<real>(ir->opts.ref_t[0], 0.0); /* using 'System' temperature */
1359 real kT = BOLTZ * reft;
1361 for (int j = 0; j < nh; j++)
1363 double iQinv = MassQ->QPinv[i * nh + j];
1366 energy += 0.5 * gmx::square(state->nhpres_vxi[i * nh + j]) / iQinv;
1367 /* contribution from the thermal variable of the NH chain */
1368 energy += state->nhpres_xi[i * nh + j] * kT;
1372 fprintf(debug, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i, j,
1373 state->nhpres_vxi[i * nh + j], state->nhpres_xi[i * nh + j]);
1381 /* Returns the energy accumulated by the V-rescale or Berendsen thermostat */
1382 static real energyVrescale(const t_inputrec* ir, const t_state* state)
1385 for (int i = 0; i < ir->opts.ngtc; i++)
1387 energy += state->therm_integral[i];
1393 real NPT_energy(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1397 if (ir->epc != epcNO)
1399 /* Compute the contribution of the pressure to the conserved quantity*/
1401 real vol = det(state->box);
1405 case epcPARRINELLORAHMAN:
1407 /* contribution from the pressure momenta */
1409 calcParrinelloRahmanInvMass(ir, state->box, invMass);
1410 for (int d = 0; d < DIM; d++)
1412 for (int n = 0; n <= d; n++)
1414 if (invMass[d][n] > 0)
1416 energyNPT += 0.5 * gmx::square(state->boxv[d][n]) / (invMass[d][n] * PRESFAC);
1421 /* Contribution from the PV term.
1422 * Not that with non-zero off-diagonal reference pressures,
1423 * i.e. applied shear stresses, there are additional terms.
1424 * We don't support this here, since that requires keeping
1425 * track of unwrapped box diagonal elements. This case is
1426 * excluded in integratorHasConservedEnergyQuantity().
1428 energyNPT += vol * trace(ir->ref_p) / (DIM * PRESFAC);
1432 /* contribution from the pressure momenta */
1433 energyNPT += 0.5 * gmx::square(state->veta) / MassQ->Winv;
1435 /* contribution from the PV term */
1436 energyNPT += vol * trace(ir->ref_p) / (DIM * PRESFAC);
1438 if (ir->epc == epcMTTK)
1440 /* contribution from the MTTK chain */
1441 energyNPT += energyPressureMTTK(ir, state, MassQ);
1444 case epcBERENDSEN: energyNPT += state->baros_integral; break;
1448 "Conserved energy quantity for pressure coupling is not handled. A case "
1449 "should be added with either the conserved quantity added or nothing added "
1450 "and an exclusion added to integratorHasConservedEnergyQuantity().");
1458 case etcBERENDSEN: energyNPT += energyVrescale(ir, state); break;
1459 case etcNOSEHOOVER: energyNPT += energyNoseHoover(ir, state, MassQ); break;
1461 case etcANDERSENMASSIVE:
1462 // Not supported, excluded in integratorHasConservedEnergyQuantity()
1467 "Conserved energy quantity for temperature coupling is not handled. A case "
1468 "should be added with either the conserved quantity added or nothing added and "
1469 "an exclusion added to integratorHasConservedEnergyQuantity().");
1476 static real vrescale_sumnoises(real nn, gmx::ThreeFry2x64<>* rng, gmx::NormalDistribution<real>* normalDist)
1479 * Returns the sum of nn independent gaussian noises squared
1480 * (i.e. equivalent to summing the square of the return values
1481 * of nn calls to a normal distribution).
1483 const real ndeg_tol = 0.0001;
1485 gmx::GammaDistribution<real> gammaDist(0.5 * nn, 1.0);
1487 if (nn < 2 + ndeg_tol)
1492 nn_int = gmx::roundToInt(nn);
1494 if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1497 "The v-rescale thermostat was called with a group with #DOF=%f, but for "
1498 "#DOF<3 only integer #DOF are supported",
1503 for (i = 0; i < nn_int; i++)
1505 gauss = (*normalDist)(*rng);
1511 /* Use a gamma distribution for any real nn > 2 */
1512 r = 2.0 * gammaDist(*rng);
1518 real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut, int64_t step, int64_t seed)
1521 * Generates a new value for the kinetic energy,
1522 * according to Bussi et al JCP (2007), Eq. (A7)
1523 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1524 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1525 * ndeg: number of degrees of freedom of the atoms to be thermalized
1526 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1528 real factor, rr, ekin_new;
1529 gmx::ThreeFry2x64<64> rng(seed, gmx::RandomDomain::Thermostat);
1530 gmx::NormalDistribution<real> normalDist;
1534 factor = exp(-1.0 / taut);
1541 rng.restart(step, 0);
1543 rr = normalDist(rng);
1547 * (sigma * (vrescale_sumnoises(ndeg - 1, &rng, &normalDist) + rr * rr) / ndeg - kk)
1548 + 2.0 * rr * std::sqrt(kk * sigma / ndeg * (1.0 - factor) * factor);
1553 void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, real dt, double therm_integral[])
1555 const t_grpopts* opts;
1557 real Ek, Ek_ref1, Ek_ref, Ek_new;
1561 for (i = 0; (i < opts->ngtc); i++)
1565 Ek = trace(ekind->tcstat[i].ekinf);
1569 Ek = trace(ekind->tcstat[i].ekinh);
1572 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1574 Ek_ref1 = 0.5 * opts->ref_t[i] * BOLTZ;
1575 Ek_ref = Ek_ref1 * opts->nrdf[i];
1577 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i], opts->tau_t[i] / dt, step, ir->ld_seed);
1579 /* Analytically Ek_new>=0, but we check for rounding errors */
1582 ekind->tcstat[i].lambda = 0.0;
1586 ekind->tcstat[i].lambda = std::sqrt(Ek_new / Ek);
1589 therm_integral[i] -= Ek_new - Ek;
1593 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n", i, Ek_ref,
1594 Ek, Ek_new, ekind->tcstat[i].lambda);
1599 ekind->tcstat[i].lambda = 1.0;
1604 void rescale_velocities(const gmx_ekindata_t* ekind, const t_mdatoms* mdatoms, int start, int end, rvec v[])
1606 unsigned short *cACC, *cTC;
1613 gmx::ArrayRef<const t_grp_tcstat> tcstat = ekind->tcstat;
1617 gmx::ArrayRef<const t_grp_acc> gstat = ekind->grpstat;
1618 cACC = mdatoms->cACC;
1622 for (n = start; n < end; n++)
1632 /* Only scale the velocity component relative to the COM velocity */
1633 rvec_sub(v[n], gstat[ga].u, vrel);
1634 lg = tcstat[gt].lambda;
1635 for (d = 0; d < DIM; d++)
1637 v[n][d] = gstat[ga].u[d] + lg * vrel[d];
1644 for (n = start; n < end; n++)
1650 lg = tcstat[gt].lambda;
1651 for (d = 0; d < DIM; d++)
1659 //! Check whether we're doing simulated annealing
1660 bool doSimulatedAnnealing(const t_inputrec* ir)
1662 for (int i = 0; i < ir->opts.ngtc; i++)
1664 /* set bSimAnn if any group is being annealed */
1665 if (ir->opts.annealing[i] != eannNO)
1673 // TODO If we keep simulated annealing, make a proper module that
1674 // does not rely on changing inputrec.
1675 bool initSimulatedAnnealing(t_inputrec* ir, gmx::Update* upd)
1677 bool doSimAnnealing = doSimulatedAnnealing(ir);
1680 update_annealing_target_temp(ir, ir->init_t, upd);
1682 return doSimAnnealing;
1685 /* set target temperatures if we are annealing */
1686 void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd)
1688 int i, j, n, npoints;
1689 real pert, thist = 0, x;
1691 for (i = 0; i < ir->opts.ngtc; i++)
1693 npoints = ir->opts.anneal_npoints[i];
1694 switch (ir->opts.annealing[i])
1696 case eannNO: continue;
1698 /* calculate time modulo the period */
1699 pert = ir->opts.anneal_time[i][npoints - 1];
1700 n = static_cast<int>(t / pert);
1701 thist = t - n * pert; /* modulo time */
1702 /* Make sure rounding didn't get us outside the interval */
1703 if (std::fabs(thist - pert) < GMX_REAL_EPS * 100)
1708 case eannSINGLE: thist = t; break;
1710 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",
1711 i, ir->opts.ngtc, npoints);
1713 /* We are doing annealing for this group if we got here,
1714 * and we have the (relative) time as thist.
1715 * calculate target temp */
1717 while ((j < npoints - 1) && (thist > (ir->opts.anneal_time[i][j + 1])))
1721 if (j < npoints - 1)
1723 /* Found our position between points j and j+1.
1724 * Interpolate: x is the amount from j+1, (1-x) from point j
1725 * First treat possible jumps in temperature as a special case.
1727 if ((ir->opts.anneal_time[i][j + 1] - ir->opts.anneal_time[i][j]) < GMX_REAL_EPS * 100)
1729 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][j + 1];
1733 x = ((thist - ir->opts.anneal_time[i][j])
1734 / (ir->opts.anneal_time[i][j + 1] - ir->opts.anneal_time[i][j]));
1736 x * ir->opts.anneal_temp[i][j + 1] + (1 - x) * ir->opts.anneal_temp[i][j];
1741 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][npoints - 1];
1745 update_temperature_constants(upd->sd(), ir);
1748 void pleaseCiteCouplingAlgorithms(FILE* fplog, const t_inputrec& ir)
1750 if (EI_DYNAMICS(ir.eI))
1752 if (ir.etc == etcBERENDSEN)
1754 please_cite(fplog, "Berendsen84a");
1756 if (ir.etc == etcVRESCALE)
1758 please_cite(fplog, "Bussi2007a");
1760 // TODO this is actually an integrator, not a coupling algorithm
1763 please_cite(fplog, "Goga2012");