2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
45 #include "gromacs/domdec/domdec_struct.h"
46 #include "gromacs/gmxlib/nrnb.h"
47 #include "gromacs/math/functions.h"
48 #include "gromacs/math/invertmatrix.h"
49 #include "gromacs/math/units.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/math/vecdump.h"
52 #include "gromacs/mdlib/expanded.h"
53 #include "gromacs/mdlib/gmx_omp_nthreads.h"
54 #include "gromacs/mdlib/stat.h"
55 #include "gromacs/mdlib/update.h"
56 #include "gromacs/mdtypes/commrec.h"
57 #include "gromacs/mdtypes/enerdata.h"
58 #include "gromacs/mdtypes/group.h"
59 #include "gromacs/mdtypes/inputrec.h"
60 #include "gromacs/mdtypes/md_enums.h"
61 #include "gromacs/mdtypes/mdatom.h"
62 #include "gromacs/mdtypes/state.h"
63 #include "gromacs/pbcutil/boxutilities.h"
64 #include "gromacs/pbcutil/pbc.h"
65 #include "gromacs/random/gammadistribution.h"
66 #include "gromacs/random/normaldistribution.h"
67 #include "gromacs/random/tabulatednormaldistribution.h"
68 #include "gromacs/random/threefry.h"
69 #include "gromacs/random/uniformrealdistribution.h"
70 #include "gromacs/utility/cstringutil.h"
71 #include "gromacs/utility/fatalerror.h"
72 #include "gromacs/utility/pleasecite.h"
73 #include "gromacs/utility/smalloc.h"
75 #define NTROTTERPARTS 3
77 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
79 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
80 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
82 #define MAX_SUZUKI_YOSHIDA_NUM 5
83 #define SUZUKI_YOSHIDA_NUM 5
85 static const double sy_const_1[] = { 1. };
86 static const double sy_const_3[] = { 0.828981543588751, -0.657963087177502, 0.828981543588751 };
87 static const double sy_const_5[] = { 0.2967324292201065, 0.2967324292201065, -0.186929716880426,
88 0.2967324292201065, 0.2967324292201065 };
90 static const double* sy_const[] = { nullptr, sy_const_1, nullptr, sy_const_3, nullptr, sy_const_5 };
93 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
97 {0.828981543588751,-0.657963087177502,0.828981543588751},
99 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
102 /* these integration routines are only referenced inside this file */
103 static void NHC_trotter(const t_grpopts* opts,
105 const gmx_ekindata_t* ekind,
111 const t_extmass* MassQ,
112 gmx_bool bEkinAveVel)
115 /* general routine for both barostat and thermostat nose hoover chains */
118 double Ekin, Efac, reft, kT, nd;
123 int mstepsi, mstepsj;
124 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
125 int nh = opts->nhchainlength;
128 mstepsi = mstepsj = ns;
130 /* if scalefac is NULL, we are doing the NHC of the barostat */
133 if (scalefac == nullptr)
138 for (i = 0; i < nvar; i++)
141 /* make it easier to iterate by selecting
142 out the sub-array that corresponds to this T group */
146 gmx::ArrayRef<const double> iQinv;
149 iQinv = gmx::arrayRefFromArray(&MassQ->QPinv[i * nh], nh);
150 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
151 reft = std::max<real>(0, opts->ref_t[0]);
152 Ekin = gmx::square(*veta) / MassQ->Winv;
156 iQinv = gmx::arrayRefFromArray(&MassQ->Qinv[i * nh], nh);
157 const t_grp_tcstat* tcstat = &ekind->tcstat[i];
159 reft = std::max<real>(0, opts->ref_t[i]);
162 Ekin = 2 * trace(tcstat->ekinf) * tcstat->ekinscalef_nhc;
166 Ekin = 2 * trace(tcstat->ekinh) * tcstat->ekinscaleh_nhc;
171 for (mi = 0; mi < mstepsi; mi++)
173 for (mj = 0; mj < mstepsj; mj++)
175 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
176 dt = sy_const[ns][mj] * dtfull / mstepsi;
178 /* compute the thermal forces */
179 GQ[0] = iQinv[0] * (Ekin - nd * kT);
181 for (j = 0; j < nh - 1; j++)
183 if (iQinv[j + 1] > 0)
185 /* we actually don't need to update here if we save the
186 state of the GQ, but it's easier to just recompute*/
187 GQ[j + 1] = iQinv[j + 1] * ((gmx::square(ivxi[j]) / iQinv[j]) - kT);
195 ivxi[nh - 1] += 0.25 * dt * GQ[nh - 1];
196 for (j = nh - 1; j > 0; j--)
198 Efac = exp(-0.125 * dt * ivxi[j]);
199 ivxi[j - 1] = Efac * (ivxi[j - 1] * Efac + 0.25 * dt * GQ[j - 1]);
202 Efac = exp(-0.5 * dt * ivxi[0]);
211 Ekin *= (Efac * Efac);
213 /* Issue - if the KE is an average of the last and the current temperatures, then we
214 might not be able to scale the kinetic energy directly with this factor. Might
215 take more bookkeeping -- have to think about this a bit more . . . */
217 GQ[0] = iQinv[0] * (Ekin - nd * kT);
219 /* update thermostat positions */
220 for (j = 0; j < nh; j++)
222 ixi[j] += 0.5 * dt * ivxi[j];
225 for (j = 0; j < nh - 1; j++)
227 Efac = exp(-0.125 * dt * ivxi[j + 1]);
228 ivxi[j] = Efac * (ivxi[j] * Efac + 0.25 * dt * GQ[j]);
229 if (iQinv[j + 1] > 0)
231 GQ[j + 1] = iQinv[j + 1] * ((gmx::square(ivxi[j]) / iQinv[j]) - kT);
238 ivxi[nh - 1] += 0.25 * dt * GQ[nh - 1];
245 static void boxv_trotter(const t_inputrec* ir,
249 const gmx_ekindata_t* ekind,
252 const t_extmass* MassQ)
259 tensor ekinmod, localpres;
261 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
262 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
265 if (ir->epct == epctSEMIISOTROPIC)
274 /* eta is in pure units. veta is in units of ps^-1. GW is in
275 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
276 taken to use only RATIOS of eta in updating the volume. */
278 /* we take the partial pressure tensors, modify the
279 kinetic energy tensor, and recovert to pressure */
281 if (ir->opts.nrdf[0] == 0)
283 gmx_fatal(FARGS, "Barostat is coupled to a T-group with no degrees of freedom\n");
285 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
286 alpha = 1.0 + DIM / (static_cast<double>(ir->opts.nrdf[0]));
287 alpha *= ekind->tcstat[0].ekinscalef_nhc;
288 msmul(ekind->ekin, alpha, ekinmod);
289 /* for now, we use Elr = 0, because if you want to get it right, you
290 really should be using PME. Maybe print a warning? */
292 pscal = calc_pres(ir->ePBC, nwall, box, ekinmod, vir, localpres) + pcorr;
295 GW = (vol * (MassQ->Winv / PRESFAC)) * (DIM * pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
297 *veta += 0.5 * dt * GW;
301 * This file implements temperature and pressure coupling algorithms:
302 * For now only the Weak coupling and the modified weak coupling.
304 * Furthermore computation of pressure and temperature is done here
308 real calc_pres(int ePBC, int nwall, const matrix box, const tensor ekin, const tensor vir, tensor pres)
313 if (ePBC == epbcNONE || (ePBC == epbcXY && nwall != 2))
319 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
320 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
324 fac = PRESFAC * 2.0 / det(box);
325 for (n = 0; (n < DIM); n++)
327 for (m = 0; (m < DIM); m++)
329 pres[n][m] = (ekin[n][m] - vir[n][m]) * fac;
335 pr_rvecs(debug, 0, "PC: pres", pres, DIM);
336 pr_rvecs(debug, 0, "PC: ekin", ekin, DIM);
337 pr_rvecs(debug, 0, "PC: vir ", vir, DIM);
338 pr_rvecs(debug, 0, "PC: box ", box, DIM);
341 return trace(pres) / DIM;
344 real calc_temp(real ekin, real nrdf)
348 return (2.0 * ekin) / (nrdf * BOLTZ);
356 /*! \brief Sets 1/mass for Parrinello-Rahman in wInv; NOTE: PRESFAC is not included, so not in GROMACS units! */
357 static void calcParrinelloRahmanInvMass(const t_inputrec* ir, const matrix box, tensor wInv)
361 /* TODO: See if we can make the mass independent of the box size */
362 maxBoxLength = std::max(box[XX][XX], box[YY][YY]);
363 maxBoxLength = std::max(maxBoxLength, box[ZZ][ZZ]);
365 for (int d = 0; d < DIM; d++)
367 for (int n = 0; n < DIM; n++)
369 wInv[d][n] = (4 * M_PI * M_PI * ir->compress[d][n]) / (3 * ir->tau_p * ir->tau_p * maxBoxLength);
374 void parrinellorahman_pcoupl(FILE* fplog,
376 const t_inputrec* ir,
386 /* This doesn't do any coordinate updating. It just
387 * integrates the box vector equations from the calculated
388 * acceleration due to pressure difference. We also compute
389 * the tensor M which is used in update to couple the particle
390 * coordinates to the box vectors.
392 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
395 * M_nk = (h') * (h' * h + h' h) * h
397 * with the dots denoting time derivatives and h is the transformation from
398 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
399 * This also goes for the pressure and M tensors - they are transposed relative
400 * to ours. Our equation thus becomes:
403 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
405 * where b is the gromacs box matrix.
406 * Our box accelerations are given by
408 * b = vol/W inv(box') * (P-ref_P) (=h')
411 real vol = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
412 real atot, arel, change;
413 tensor invbox, pdiff, t1, t2;
415 gmx::invertBoxMatrix(box, invbox);
419 /* Note that PRESFAC does not occur here.
420 * The pressure and compressibility always occur as a product,
421 * therefore the pressure unit drops out.
424 calcParrinelloRahmanInvMass(ir, box, winv);
426 m_sub(pres, ir->ref_p, pdiff);
428 if (ir->epct == epctSURFACETENSION)
430 /* Unlike Berendsen coupling it might not be trivial to include a z
431 * pressure correction here? On the other hand we don't scale the
432 * box momentarily, but change accelerations, so it might not be crucial.
434 real xy_pressure = 0.5 * (pres[XX][XX] + pres[YY][YY]);
435 for (int d = 0; d < ZZ; d++)
437 pdiff[d][d] = (xy_pressure - (pres[ZZ][ZZ] - ir->ref_p[d][d] / box[d][d]));
441 tmmul(invbox, pdiff, t1);
442 /* Move the off-diagonal elements of the 'force' to one side to ensure
443 * that we obey the box constraints.
445 for (int d = 0; d < DIM; d++)
447 for (int n = 0; n < d; n++)
449 t1[d][n] += t1[n][d];
456 case epctANISOTROPIC:
457 for (int d = 0; d < DIM; d++)
459 for (int n = 0; n <= d; n++)
461 t1[d][n] *= winv[d][n] * vol;
466 /* calculate total volume acceleration */
467 atot = box[XX][XX] * box[YY][YY] * t1[ZZ][ZZ] + box[XX][XX] * t1[YY][YY] * box[ZZ][ZZ]
468 + t1[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
469 arel = atot / (3 * vol);
470 /* set all RELATIVE box accelerations equal, and maintain total V
472 for (int d = 0; d < DIM; d++)
474 for (int n = 0; n <= d; n++)
476 t1[d][n] = winv[0][0] * vol * arel * box[d][n];
480 case epctSEMIISOTROPIC:
481 case epctSURFACETENSION:
482 /* Note the correction to pdiff above for surftens. coupling */
484 /* calculate total XY volume acceleration */
485 atot = box[XX][XX] * t1[YY][YY] + t1[XX][XX] * box[YY][YY];
486 arel = atot / (2 * box[XX][XX] * box[YY][YY]);
487 /* set RELATIVE XY box accelerations equal, and maintain total V
488 * change speed. Dont change the third box vector accelerations */
489 for (int d = 0; d < ZZ; d++)
491 for (int n = 0; n <= d; n++)
493 t1[d][n] = winv[d][n] * vol * arel * box[d][n];
496 for (int n = 0; n < DIM; n++)
498 t1[ZZ][n] *= winv[ZZ][n] * vol;
503 "Parrinello-Rahman pressure coupling type %s "
504 "not supported yet\n",
505 EPCOUPLTYPETYPE(ir->epct));
509 for (int d = 0; d < DIM; d++)
511 for (int n = 0; n <= d; n++)
513 boxv[d][n] += dt * t1[d][n];
515 /* We do NOT update the box vectors themselves here, since
516 * we need them for shifting later. It is instead done last
517 * in the update() routine.
520 /* Calculate the change relative to diagonal elements-
521 since it's perfectly ok for the off-diagonal ones to
522 be zero it doesn't make sense to check the change relative
526 change = std::fabs(dt * boxv[d][n] / box[d][d]);
528 if (change > maxchange)
535 if (maxchange > 0.01 && fplog)
539 "\nStep %s Warning: Pressure scaling more than 1%%. "
540 "This may mean your system\n is not yet equilibrated. "
541 "Use of Parrinello-Rahman pressure coupling during\n"
542 "equilibration can lead to simulation instability, "
543 "and is discouraged.\n",
544 gmx_step_str(step, buf));
548 preserve_box_shape(ir, box_rel, boxv);
550 mtmul(boxv, box, t1); /* t1=boxv * b' */
551 mmul(invbox, t1, t2);
552 mtmul(t2, invbox, M);
554 /* Determine the scaling matrix mu for the coordinates */
555 for (int d = 0; d < DIM; d++)
557 for (int n = 0; n <= d; n++)
559 t1[d][n] = box[d][n] + dt * boxv[d][n];
562 preserve_box_shape(ir, box_rel, t1);
563 /* t1 is the box at t+dt, determine mu as the relative change */
564 mmul_ur0(invbox, t1, mu);
567 void berendsen_pcoupl(FILE* fplog,
569 const t_inputrec* ir,
573 const matrix force_vir,
574 const matrix constraint_vir,
576 double* baros_integral)
579 real scalar_pressure, xy_pressure, p_corr_z;
583 * Calculate the scaling matrix mu
587 for (d = 0; d < DIM; d++)
589 scalar_pressure += pres[d][d] / DIM;
592 xy_pressure += pres[d][d] / (DIM - 1);
595 /* Pressure is now in bar, everywhere. */
596 #define factor(d, m) (ir->compress[d][m] * dt / ir->tau_p)
598 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
599 * necessary for triclinic scaling
605 for (d = 0; d < DIM; d++)
607 mu[d][d] = 1.0 - factor(d, d) * (ir->ref_p[d][d] - scalar_pressure) / DIM;
610 case epctSEMIISOTROPIC:
611 for (d = 0; d < ZZ; d++)
613 mu[d][d] = 1.0 - factor(d, d) * (ir->ref_p[d][d] - xy_pressure) / DIM;
615 mu[ZZ][ZZ] = 1.0 - factor(ZZ, ZZ) * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM;
617 case epctANISOTROPIC:
618 for (d = 0; d < DIM; d++)
620 for (n = 0; n < DIM; n++)
622 mu[d][n] = (d == n ? 1.0 : 0.0) - factor(d, n) * (ir->ref_p[d][n] - pres[d][n]) / DIM;
626 case epctSURFACETENSION:
627 /* ir->ref_p[0/1] is the reference surface-tension times *
628 * the number of surfaces */
629 if (ir->compress[ZZ][ZZ] != 0.0F)
631 p_corr_z = dt / ir->tau_p * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
635 /* when the compressibity is zero, set the pressure correction *
636 * in the z-direction to zero to get the correct surface tension */
639 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ] * p_corr_z;
640 for (d = 0; d < DIM - 1; d++)
644 * (ir->ref_p[d][d] / (mu[ZZ][ZZ] * box[ZZ][ZZ])
645 - (pres[ZZ][ZZ] + p_corr_z - xy_pressure))
650 gmx_fatal(FARGS, "Berendsen pressure coupling type %s not supported yet\n",
651 EPCOUPLTYPETYPE(ir->epct));
653 /* To fullfill the orientation restrictions on triclinic boxes
654 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
655 * the other elements of mu to first order.
657 mu[YY][XX] += mu[XX][YY];
658 mu[ZZ][XX] += mu[XX][ZZ];
659 mu[ZZ][YY] += mu[YY][ZZ];
664 /* Keep track of the work the barostat applies on the system.
665 * Without constraints force_vir tells us how Epot changes when scaling.
666 * With constraints constraint_vir gives us the constraint contribution
667 * to both Epot and Ekin. Although we are not scaling velocities, scaling
668 * the coordinates leads to scaling of distances involved in constraints.
669 * This in turn changes the angular momentum (even if the constrained
670 * distances are corrected at the next step). The kinetic component
671 * of the constraint virial captures the angular momentum change.
673 for (int d = 0; d < DIM; d++)
675 for (int n = 0; n <= d; n++)
678 2 * (mu[d][n] - (n == d ? 1 : 0)) * (force_vir[d][n] + constraint_vir[d][n]);
684 pr_rvecs(debug, 0, "PC: pres ", pres, 3);
685 pr_rvecs(debug, 0, "PC: mu ", mu, 3);
688 if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 || mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01
689 || mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
693 "\nStep %s Warning: pressure scaling more than 1%%, "
695 gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
698 fprintf(fplog, "%s", buf);
700 fprintf(stderr, "%s", buf);
704 void berendsen_pscale(const t_inputrec* ir,
711 const unsigned short cFREEZE[],
713 const bool scaleCoordinates)
715 ivec* nFreeze = ir->opts.nFreeze;
717 int nthreads gmx_unused;
719 #ifndef __clang_analyzer__
720 nthreads = gmx_omp_nthreads_get(emntUpdate);
723 /* Scale the positions */
724 if (scaleCoordinates)
726 #pragma omp parallel for num_threads(nthreads) schedule(static)
727 for (int n = start; n < start + nr_atoms; n++)
729 // Trivial OpenMP region that does not throw
732 if (cFREEZE == nullptr)
743 x[n][XX] = mu[XX][XX] * x[n][XX] + mu[YY][XX] * x[n][YY] + mu[ZZ][XX] * x[n][ZZ];
747 x[n][YY] = mu[YY][YY] * x[n][YY] + mu[ZZ][YY] * x[n][ZZ];
751 x[n][ZZ] = mu[ZZ][ZZ] * x[n][ZZ];
755 /* compute final boxlengths */
756 for (d = 0; d < DIM; d++)
758 box[d][XX] = mu[XX][XX] * box[d][XX] + mu[YY][XX] * box[d][YY] + mu[ZZ][XX] * box[d][ZZ];
759 box[d][YY] = mu[YY][YY] * box[d][YY] + mu[ZZ][YY] * box[d][ZZ];
760 box[d][ZZ] = mu[ZZ][ZZ] * box[d][ZZ];
763 preserve_box_shape(ir, box_rel, box);
765 /* (un)shifting should NOT be done after this,
766 * since the box vectors might have changed
768 inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
771 void berendsen_tcoupl(const t_inputrec* ir, gmx_ekindata_t* ekind, real dt, std::vector<double>& therm_integral)
773 const t_grpopts* opts = &ir->opts;
775 for (int i = 0; (i < opts->ngtc); i++)
781 Ek = trace(ekind->tcstat[i].ekinf);
782 T = ekind->tcstat[i].T;
786 Ek = trace(ekind->tcstat[i].ekinh);
787 T = ekind->tcstat[i].Th;
790 if ((opts->tau_t[i] > 0) && (T > 0.0))
792 real reft = std::max<real>(0, opts->ref_t[i]);
793 real lll = std::sqrt(1.0 + (dt / opts->tau_t[i]) * (reft / T - 1.0));
794 ekind->tcstat[i].lambda = std::max<real>(std::min<real>(lll, 1.25), 0.8);
798 ekind->tcstat[i].lambda = 1.0;
801 /* Keep track of the amount of energy we are adding to the system */
802 therm_integral[i] -= (gmx::square(ekind->tcstat[i].lambda) - 1) * Ek;
806 fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n", i, T, ekind->tcstat[i].lambda);
811 void andersen_tcoupl(const t_inputrec* ir,
815 gmx::ArrayRef<gmx::RVec> v,
817 const std::vector<bool>& randomize,
818 gmx::ArrayRef<const real> boltzfac)
820 const int* gatindex = (DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
823 gmx::ThreeFry2x64<0> rng(ir->andersen_seed, gmx::RandomDomain::Thermostat);
824 gmx::UniformRealDistribution<real> uniformDist;
825 gmx::TabulatedNormalDistribution<real, 14> normalDist;
827 /* randomize the velocities of the selected particles */
829 for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */
831 int ng = gatindex ? gatindex[i] : i;
834 rng.restart(step, ng);
838 gc = md->cTC[i]; /* assign the atom to a temperature group if there are more than one */
842 if (ir->etc == etcANDERSENMASSIVE)
844 /* Randomize particle always */
849 /* Randomize particle probabilistically */
851 bRandomize = uniformDist(rng) < rate;
858 scal = std::sqrt(boltzfac[gc] * md->invmass[i]);
862 for (d = 0; d < DIM; d++)
864 v[i][d] = scal * normalDist(rng);
872 void nosehoover_tcoupl(const t_grpopts* opts,
873 const gmx_ekindata_t* ekind,
877 const t_extmass* MassQ)
882 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
884 for (i = 0; (i < opts->ngtc); i++)
886 reft = std::max<real>(0, opts->ref_t[i]);
888 vxi[i] += dt * MassQ->Qinv[i] * (ekind->tcstat[i].Th - reft);
889 xi[i] += dt * (oldvxi + vxi[i]) * 0.5;
893 void trotter_update(const t_inputrec* ir,
895 gmx_ekindata_t* ekind,
896 const gmx_enerdata_t* enerd,
900 const t_extmass* MassQ,
901 gmx::ArrayRef<std::vector<int>> trotter_seqlist,
905 int n, i, d, ngtc, gc = 0, t;
906 t_grp_tcstat* tcstat;
907 const t_grpopts* opts;
910 double * scalefac, dtc;
911 rvec sumv = { 0, 0, 0 };
914 if (trotter_seqno <= ettTSEQ2)
916 step_eff = step - 1; /* the velocity verlet calls are actually out of order -- the first
917 half step is actually the last half step from the previous step.
918 Thus the first half step actually corresponds to the n-1 step*/
925 bCouple = (ir->nsttcouple == 1 || do_per_step(step_eff + ir->nsttcouple, ir->nsttcouple));
927 const gmx::ArrayRef<const int> trotter_seq = trotter_seqlist[trotter_seqno];
929 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
933 dtc = ir->nsttcouple * ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
934 opts = &(ir->opts); /* just for ease of referencing */
937 snew(scalefac, opts->ngtc);
938 for (i = 0; i < ngtc; i++)
942 /* execute the series of trotter updates specified in the trotterpart array */
944 for (i = 0; i < NTROTTERPARTS; i++)
946 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
947 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2)
948 || (trotter_seq[i] == etrtNHC2))
957 auto v = makeArrayRef(state->v);
958 switch (trotter_seq[i])
962 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir,
963 enerd->term[F_PDISPCORR], MassQ);
967 NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi.data(),
968 state->nhpres_vxi.data(), nullptr, &(state->veta), MassQ, FALSE);
972 NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi.data(),
973 state->nosehoover_vxi.data(), scalefac, nullptr, MassQ, (ir->eI == eiVV));
974 /* need to rescale the kinetic energies and velocities here. Could
975 scale the velocities later, but we need them scaled in order to
976 produce the correct outputs, so we'll scale them here. */
978 for (t = 0; t < ngtc; t++)
980 tcstat = &ekind->tcstat[t];
981 tcstat->vscale_nhc = scalefac[t];
982 tcstat->ekinscaleh_nhc *= (scalefac[t] * scalefac[t]);
983 tcstat->ekinscalef_nhc *= (scalefac[t] * scalefac[t]);
985 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
986 /* but do we actually need the total? */
988 /* modify the velocities as well */
989 for (n = 0; n < md->homenr; n++)
991 if (md->cTC) /* does this conditional need to be here? is this always true?*/
995 for (d = 0; d < DIM; d++)
997 v[n][d] *= scalefac[gc];
1002 for (d = 0; d < DIM; d++)
1004 sumv[d] += (v[n][d]) / md->invmass[n];
1012 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
1017 extern void init_npt_masses(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool bInit)
1019 int n, i, j, d, ngtc, nh;
1020 const t_grpopts* opts;
1021 real reft, kT, ndj, nd;
1023 opts = &(ir->opts); /* just for ease of referencing */
1024 ngtc = ir->opts.ngtc;
1025 nh = state->nhchainlength;
1031 MassQ->Qinv.resize(ngtc);
1033 for (i = 0; (i < ngtc); i++)
1035 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1037 MassQ->Qinv[i] = 1.0 / (gmx::square(opts->tau_t[i] / M_2PI) * opts->ref_t[i]);
1041 MassQ->Qinv[i] = 0.0;
1045 else if (EI_VV(ir->eI))
1047 /* Set pressure variables */
1051 if (state->vol0 == 0)
1053 state->vol0 = det(state->box);
1054 /* because we start by defining a fixed
1055 compressibility, we need the volume at this
1056 compressibility to solve the problem. */
1060 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1061 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1062 MassQ->Winv = (PRESFAC * trace(ir->compress) * BOLTZ * opts->ref_t[0])
1063 / (DIM * state->vol0 * gmx::square(ir->tau_p / M_2PI));
1064 /* An alternate mass definition, from Tuckerman et al. */
1065 /* MassQ->Winv = 1.0/(gmx::square(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1066 for (d = 0; d < DIM; d++)
1068 for (n = 0; n < DIM; n++)
1070 MassQ->Winvm[d][n] =
1071 PRESFAC * ir->compress[d][n] / (state->vol0 * gmx::square(ir->tau_p / M_2PI));
1072 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1073 before using MTTK for anisotropic states.*/
1076 /* Allocate space for thermostat variables */
1079 MassQ->Qinv.resize(ngtc * nh);
1082 /* now, set temperature variables */
1083 for (i = 0; i < ngtc; i++)
1085 if (opts->tau_t[i] > 0 && opts->ref_t[i] > 0 && opts->nrdf[i] > 0)
1087 reft = std::max<real>(0, opts->ref_t[i]);
1090 for (j = 0; j < nh; j++)
1100 MassQ->Qinv[i * nh + j] = 1.0 / (gmx::square(opts->tau_t[i] / M_2PI) * ndj * kT);
1105 for (j = 0; j < nh; j++)
1107 MassQ->Qinv[i * nh + j] = 0.0;
1114 std::array<std::vector<int>, ettTSEQMAX>
1115 init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool bTrotter)
1117 int i, j, nnhpres, nh;
1118 const t_grpopts* opts;
1119 real bmass, qmass, reft, kT;
1121 opts = &(ir->opts); /* just for ease of referencing */
1122 nnhpres = state->nnhpres;
1123 nh = state->nhchainlength;
1125 if (EI_VV(ir->eI) && (ir->epc == epcMTTK) && (ir->etc != etcNOSEHOOVER))
1128 "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
1131 init_npt_masses(ir, state, MassQ, TRUE);
1133 /* first, initialize clear all the trotter calls */
1134 std::array<std::vector<int>, ettTSEQMAX> trotter_seq;
1135 for (i = 0; i < ettTSEQMAX; i++)
1137 trotter_seq[i].resize(NTROTTERPARTS, etrtNONE);
1138 trotter_seq[i][0] = etrtSKIPALL;
1143 /* no trotter calls, so we never use the values in the array.
1144 * We access them (so we need to define them, but ignore
1150 /* compute the kinetic energy by using the half step velocities or
1151 * the kinetic energies, depending on the order of the trotter calls */
1155 if (inputrecNptTrotter(ir))
1157 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1158 We start with the initial one. */
1159 /* first, a round that estimates veta. */
1160 trotter_seq[0][0] = etrtBAROV;
1162 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1164 /* The first half trotter update */
1165 trotter_seq[2][0] = etrtBAROV;
1166 trotter_seq[2][1] = etrtNHC;
1167 trotter_seq[2][2] = etrtBARONHC;
1169 /* The second half trotter update */
1170 trotter_seq[3][0] = etrtBARONHC;
1171 trotter_seq[3][1] = etrtNHC;
1172 trotter_seq[3][2] = etrtBAROV;
1174 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1176 else if (inputrecNvtTrotter(ir))
1178 /* This is the easy version - there are only two calls, both the same.
1179 Otherwise, even easier -- no calls */
1180 trotter_seq[2][0] = etrtNHC;
1181 trotter_seq[3][0] = etrtNHC;
1183 else if (inputrecNphTrotter(ir))
1185 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1186 We start with the initial one. */
1187 /* first, a round that estimates veta. */
1188 trotter_seq[0][0] = etrtBAROV;
1190 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1192 /* The first half trotter update */
1193 trotter_seq[2][0] = etrtBAROV;
1194 trotter_seq[2][1] = etrtBARONHC;
1196 /* The second half trotter update */
1197 trotter_seq[3][0] = etrtBARONHC;
1198 trotter_seq[3][1] = etrtBAROV;
1200 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1203 else if (ir->eI == eiVVAK)
1205 if (inputrecNptTrotter(ir))
1207 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1208 We start with the initial one. */
1209 /* first, a round that estimates veta. */
1210 trotter_seq[0][0] = etrtBAROV;
1212 /* The first half trotter update, part 1 -- double update, because it commutes */
1213 trotter_seq[1][0] = etrtNHC;
1215 /* The first half trotter update, part 2 */
1216 trotter_seq[2][0] = etrtBAROV;
1217 trotter_seq[2][1] = etrtBARONHC;
1219 /* The second half trotter update, part 1 */
1220 trotter_seq[3][0] = etrtBARONHC;
1221 trotter_seq[3][1] = etrtBAROV;
1223 /* The second half trotter update */
1224 trotter_seq[4][0] = etrtNHC;
1226 else if (inputrecNvtTrotter(ir))
1228 /* This is the easy version - there is only one call, both the same.
1229 Otherwise, even easier -- no calls */
1230 trotter_seq[1][0] = etrtNHC;
1231 trotter_seq[4][0] = etrtNHC;
1233 else if (inputrecNphTrotter(ir))
1235 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1236 We start with the initial one. */
1237 /* first, a round that estimates veta. */
1238 trotter_seq[0][0] = etrtBAROV;
1240 /* The first half trotter update, part 1 -- leave zero */
1241 trotter_seq[1][0] = etrtNHC;
1243 /* The first half trotter update, part 2 */
1244 trotter_seq[2][0] = etrtBAROV;
1245 trotter_seq[2][1] = etrtBARONHC;
1247 /* The second half trotter update, part 1 */
1248 trotter_seq[3][0] = etrtBARONHC;
1249 trotter_seq[3][1] = etrtBAROV;
1251 /* The second half trotter update -- blank for now */
1258 default: bmass = DIM * DIM; /* recommended mass parameters for isotropic barostat */
1261 MassQ->QPinv.resize(nnhpres * opts->nhchainlength);
1263 /* barostat temperature */
1264 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1266 reft = std::max<real>(0, opts->ref_t[0]);
1268 for (i = 0; i < nnhpres; i++)
1270 for (j = 0; j < nh; j++)
1280 MassQ->QPinv[i * opts->nhchainlength + j] =
1281 1.0 / (gmx::square(opts->tau_t[0] / M_2PI) * qmass * kT);
1287 for (i = 0; i < nnhpres; i++)
1289 for (j = 0; j < nh; j++)
1291 MassQ->QPinv[i * nh + j] = 0.0;
1298 static real energyNoseHoover(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1302 int nh = state->nhchainlength;
1304 for (int i = 0; i < ir->opts.ngtc; i++)
1306 const double* ixi = &state->nosehoover_xi[i * nh];
1307 const double* ivxi = &state->nosehoover_vxi[i * nh];
1308 const double* iQinv = &(MassQ->Qinv[i * nh]);
1310 int nd = static_cast<int>(ir->opts.nrdf[i]);
1311 real reft = std::max<real>(ir->opts.ref_t[i], 0);
1312 real kT = BOLTZ * reft;
1316 if (inputrecNvtTrotter(ir))
1318 /* contribution from the thermal momenta of the NH chain */
1319 for (int j = 0; j < nh; j++)
1323 energy += 0.5 * gmx::square(ivxi[j]) / iQinv[j];
1324 /* contribution from the thermal variable of the NH chain */
1334 energy += ndj * ixi[j] * kT;
1338 else /* Other non Trotter temperature NH control -- no chains yet. */
1340 energy += 0.5 * BOLTZ * nd * gmx::square(ivxi[0]) / iQinv[0];
1341 energy += nd * ixi[0] * kT;
1349 /* Returns the energy from the barostat thermostat chain */
1350 static real energyPressureMTTK(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1354 int nh = state->nhchainlength;
1356 for (int i = 0; i < state->nnhpres; i++)
1358 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1359 real reft = std::max<real>(ir->opts.ref_t[0], 0.0); /* using 'System' temperature */
1360 real kT = BOLTZ * reft;
1362 for (int j = 0; j < nh; j++)
1364 double iQinv = MassQ->QPinv[i * nh + j];
1367 energy += 0.5 * gmx::square(state->nhpres_vxi[i * nh + j] / iQinv);
1368 /* contribution from the thermal variable of the NH chain */
1369 energy += state->nhpres_xi[i * nh + j] * kT;
1373 fprintf(debug, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i, j,
1374 state->nhpres_vxi[i * nh + j], state->nhpres_xi[i * nh + j]);
1382 /* Returns the energy accumulated by the V-rescale or Berendsen thermostat */
1383 static real energyVrescale(const t_inputrec* ir, const t_state* state)
1386 for (int i = 0; i < ir->opts.ngtc; i++)
1388 energy += state->therm_integral[i];
1394 real NPT_energy(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1398 if (ir->epc != epcNO)
1400 /* Compute the contribution of the pressure to the conserved quantity*/
1402 real vol = det(state->box);
1406 case epcPARRINELLORAHMAN:
1408 /* contribution from the pressure momenta */
1410 calcParrinelloRahmanInvMass(ir, state->box, invMass);
1411 for (int d = 0; d < DIM; d++)
1413 for (int n = 0; n <= d; n++)
1415 if (invMass[d][n] > 0)
1417 energyNPT += 0.5 * gmx::square(state->boxv[d][n]) / (invMass[d][n] * PRESFAC);
1422 /* Contribution from the PV term.
1423 * Not that with non-zero off-diagonal reference pressures,
1424 * i.e. applied shear stresses, there are additional terms.
1425 * We don't support this here, since that requires keeping
1426 * track of unwrapped box diagonal elements. This case is
1427 * excluded in integratorHasConservedEnergyQuantity().
1429 energyNPT += vol * trace(ir->ref_p) / (DIM * PRESFAC);
1433 /* contribution from the pressure momenta */
1434 energyNPT += 0.5 * gmx::square(state->veta) / MassQ->Winv;
1436 /* contribution from the PV term */
1437 energyNPT += vol * trace(ir->ref_p) / (DIM * PRESFAC);
1439 if (ir->epc == epcMTTK)
1441 /* contribution from the MTTK chain */
1442 energyNPT += energyPressureMTTK(ir, state, MassQ);
1445 case epcBERENDSEN: energyNPT += state->baros_integral; break;
1449 "Conserved energy quantity for pressure coupling is not handled. A case "
1450 "should be added with either the conserved quantity added or nothing added "
1451 "and an exclusion added to integratorHasConservedEnergyQuantity().");
1459 case etcBERENDSEN: energyNPT += energyVrescale(ir, state); break;
1460 case etcNOSEHOOVER: energyNPT += energyNoseHoover(ir, state, MassQ); break;
1462 case etcANDERSENMASSIVE:
1463 // Not supported, excluded in integratorHasConservedEnergyQuantity()
1468 "Conserved energy quantity for temperature coupling is not handled. A case "
1469 "should be added with either the conserved quantity added or nothing added and "
1470 "an exclusion added to integratorHasConservedEnergyQuantity().");
1477 static real vrescale_sumnoises(real nn, gmx::ThreeFry2x64<>* rng, gmx::NormalDistribution<real>* normalDist)
1480 * Returns the sum of nn independent gaussian noises squared
1481 * (i.e. equivalent to summing the square of the return values
1482 * of nn calls to a normal distribution).
1484 const real ndeg_tol = 0.0001;
1486 gmx::GammaDistribution<real> gammaDist(0.5 * nn, 1.0);
1488 if (nn < 2 + ndeg_tol)
1493 nn_int = gmx::roundToInt(nn);
1495 if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1498 "The v-rescale thermostat was called with a group with #DOF=%f, but for "
1499 "#DOF<3 only integer #DOF are supported",
1504 for (i = 0; i < nn_int; i++)
1506 gauss = (*normalDist)(*rng);
1512 /* Use a gamma distribution for any real nn > 2 */
1513 r = 2.0 * gammaDist(*rng);
1519 real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut, int64_t step, int64_t seed)
1522 * Generates a new value for the kinetic energy,
1523 * according to Bussi et al JCP (2007), Eq. (A7)
1524 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1525 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1526 * ndeg: number of degrees of freedom of the atoms to be thermalized
1527 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1529 real factor, rr, ekin_new;
1530 gmx::ThreeFry2x64<64> rng(seed, gmx::RandomDomain::Thermostat);
1531 gmx::NormalDistribution<real> normalDist;
1535 factor = exp(-1.0 / taut);
1542 rng.restart(step, 0);
1544 rr = normalDist(rng);
1548 * (sigma * (vrescale_sumnoises(ndeg - 1, &rng, &normalDist) + rr * rr) / ndeg - kk)
1549 + 2.0 * rr * std::sqrt(kk * sigma / ndeg * (1.0 - factor) * factor);
1554 void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, real dt, double therm_integral[])
1556 const t_grpopts* opts;
1558 real Ek, Ek_ref1, Ek_ref, Ek_new;
1562 for (i = 0; (i < opts->ngtc); i++)
1566 Ek = trace(ekind->tcstat[i].ekinf);
1570 Ek = trace(ekind->tcstat[i].ekinh);
1573 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1575 Ek_ref1 = 0.5 * opts->ref_t[i] * BOLTZ;
1576 Ek_ref = Ek_ref1 * opts->nrdf[i];
1578 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i], opts->tau_t[i] / dt, step, ir->ld_seed);
1580 /* Analytically Ek_new>=0, but we check for rounding errors */
1583 ekind->tcstat[i].lambda = 0.0;
1587 ekind->tcstat[i].lambda = std::sqrt(Ek_new / Ek);
1590 therm_integral[i] -= Ek_new - Ek;
1594 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n", i, Ek_ref,
1595 Ek, Ek_new, ekind->tcstat[i].lambda);
1600 ekind->tcstat[i].lambda = 1.0;
1605 void rescale_velocities(const gmx_ekindata_t* ekind, const t_mdatoms* mdatoms, int start, int end, rvec v[])
1607 unsigned short *cACC, *cTC;
1614 gmx::ArrayRef<const t_grp_tcstat> tcstat = ekind->tcstat;
1618 gmx::ArrayRef<const t_grp_acc> gstat = ekind->grpstat;
1619 cACC = mdatoms->cACC;
1623 for (n = start; n < end; n++)
1633 /* Only scale the velocity component relative to the COM velocity */
1634 rvec_sub(v[n], gstat[ga].u, vrel);
1635 lg = tcstat[gt].lambda;
1636 for (d = 0; d < DIM; d++)
1638 v[n][d] = gstat[ga].u[d] + lg * vrel[d];
1645 for (n = start; n < end; n++)
1651 lg = tcstat[gt].lambda;
1652 for (d = 0; d < DIM; d++)
1660 //! Check whether we're doing simulated annealing
1661 bool doSimulatedAnnealing(const t_inputrec* ir)
1663 for (int i = 0; i < ir->opts.ngtc; i++)
1665 /* set bSimAnn if any group is being annealed */
1666 if (ir->opts.annealing[i] != eannNO)
1674 // TODO If we keep simulated annealing, make a proper module that
1675 // does not rely on changing inputrec.
1676 bool initSimulatedAnnealing(t_inputrec* ir, gmx::Update* upd)
1678 bool doSimAnnealing = doSimulatedAnnealing(ir);
1681 update_annealing_target_temp(ir, ir->init_t, upd);
1683 return doSimAnnealing;
1686 /* set target temperatures if we are annealing */
1687 void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd)
1689 int i, j, n, npoints;
1690 real pert, thist = 0, x;
1692 for (i = 0; i < ir->opts.ngtc; i++)
1694 npoints = ir->opts.anneal_npoints[i];
1695 switch (ir->opts.annealing[i])
1697 case eannNO: continue;
1699 /* calculate time modulo the period */
1700 pert = ir->opts.anneal_time[i][npoints - 1];
1701 n = static_cast<int>(t / pert);
1702 thist = t - n * pert; /* modulo time */
1703 /* Make sure rounding didn't get us outside the interval */
1704 if (std::fabs(thist - pert) < GMX_REAL_EPS * 100)
1709 case eannSINGLE: thist = t; break;
1711 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",
1712 i, ir->opts.ngtc, npoints);
1714 /* We are doing annealing for this group if we got here,
1715 * and we have the (relative) time as thist.
1716 * calculate target temp */
1718 while ((j < npoints - 1) && (thist > (ir->opts.anneal_time[i][j + 1])))
1722 if (j < npoints - 1)
1724 /* Found our position between points j and j+1.
1725 * Interpolate: x is the amount from j+1, (1-x) from point j
1726 * First treat possible jumps in temperature as a special case.
1728 if ((ir->opts.anneal_time[i][j + 1] - ir->opts.anneal_time[i][j]) < GMX_REAL_EPS * 100)
1730 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][j + 1];
1734 x = ((thist - ir->opts.anneal_time[i][j])
1735 / (ir->opts.anneal_time[i][j + 1] - ir->opts.anneal_time[i][j]));
1737 x * ir->opts.anneal_temp[i][j + 1] + (1 - x) * ir->opts.anneal_temp[i][j];
1742 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][npoints - 1];
1746 update_temperature_constants(upd->sd(), ir);
1749 void pleaseCiteCouplingAlgorithms(FILE* fplog, const t_inputrec& ir)
1751 if (EI_DYNAMICS(ir.eI))
1753 if (ir.etc == etcBERENDSEN)
1755 please_cite(fplog, "Berendsen84a");
1757 if (ir.etc == etcVRESCALE)
1759 please_cite(fplog, "Bussi2007a");
1761 // TODO this is actually an integrator, not a coupling algorithm
1764 please_cite(fplog, "Goga2012");