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,2021, 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.
47 #include "gromacs/domdec/domdec_struct.h"
48 #include "gromacs/gmxlib/nrnb.h"
49 #include "gromacs/math/functions.h"
50 #include "gromacs/math/invertmatrix.h"
51 #include "gromacs/math/units.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/math/vecdump.h"
54 #include "gromacs/mdlib/boxdeformation.h"
55 #include "gromacs/mdlib/expanded.h"
56 #include "gromacs/mdlib/gmx_omp_nthreads.h"
57 #include "gromacs/mdlib/stat.h"
58 #include "gromacs/mdlib/update.h"
59 #include "gromacs/mdtypes/commrec.h"
60 #include "gromacs/mdtypes/enerdata.h"
61 #include "gromacs/mdtypes/group.h"
62 #include "gromacs/mdtypes/inputrec.h"
63 #include "gromacs/mdtypes/md_enums.h"
64 #include "gromacs/mdtypes/state.h"
65 #include "gromacs/pbcutil/boxutilities.h"
66 #include "gromacs/pbcutil/pbc.h"
67 #include "gromacs/random/gammadistribution.h"
68 #include "gromacs/random/normaldistribution.h"
69 #include "gromacs/random/tabulatednormaldistribution.h"
70 #include "gromacs/random/threefry.h"
71 #include "gromacs/random/uniformrealdistribution.h"
72 #include "gromacs/utility/cstringutil.h"
73 #include "gromacs/utility/fatalerror.h"
74 #include "gromacs/utility/pleasecite.h"
75 #include "gromacs/utility/smalloc.h"
77 #define NTROTTERPARTS 3
79 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
81 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
82 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
84 #define MAX_SUZUKI_YOSHIDA_NUM 5
85 #define SUZUKI_YOSHIDA_NUM 5
87 static const double sy_const_1[] = { 1. };
88 static const double sy_const_3[] = { 0.828981543588751, -0.657963087177502, 0.828981543588751 };
89 static const double sy_const_5[] = { 0.2967324292201065,
95 static constexpr std::array<const double*, 6> sy_const = { nullptr, sy_const_1, nullptr,
96 sy_const_3, nullptr, sy_const_5 };
99 void update_tcouple(int64_t step,
100 const t_inputrec* inputrec,
102 gmx_ekindata_t* ekind,
103 const t_extmass* MassQ,
105 gmx::ArrayRef<const unsigned short> cTC)
108 // This condition was explicitly checked in previous version, but should have never been satisfied
109 GMX_ASSERT(!(EI_VV(inputrec->eI)
110 && (inputrecNvtTrotter(inputrec) || inputrecNptTrotter(inputrec)
111 || inputrecNphTrotter(inputrec))),
112 "Temperature coupling was requested with velocity verlet and trotter");
114 bool doTemperatureCoupling = false;
116 // For VV temperature coupling parameters are updated on the current
117 // step, for the others - one step before.
118 if (inputrec->etc == TemperatureCoupling::No)
120 doTemperatureCoupling = false;
122 else if (EI_VV(inputrec->eI))
124 doTemperatureCoupling = do_per_step(step, inputrec->nsttcouple);
128 doTemperatureCoupling = do_per_step(step + inputrec->nsttcouple - 1, inputrec->nsttcouple);
131 if (doTemperatureCoupling)
133 real dttc = inputrec->nsttcouple * inputrec->delta_t;
135 // TODO: berendsen_tcoupl(...), nosehoover_tcoupl(...) and vrescale_tcoupl(...) update
136 // temperature coupling parameters, which should be reflected in the name of these
138 switch (inputrec->etc)
140 case TemperatureCoupling::No:
141 case TemperatureCoupling::Andersen:
142 case TemperatureCoupling::AndersenMassive: break;
143 case TemperatureCoupling::Berendsen:
144 berendsen_tcoupl(inputrec, ekind, dttc, state->therm_integral);
146 case TemperatureCoupling::NoseHoover:
148 &(inputrec->opts), ekind, dttc, state->nosehoover_xi, state->nosehoover_vxi, MassQ);
150 case TemperatureCoupling::VRescale:
151 vrescale_tcoupl(inputrec, step, ekind, dttc, state->therm_integral);
153 default: gmx_fatal(FARGS, "Unknown temperature coupling algorithm");
155 /* rescale in place here */
156 if (EI_VV(inputrec->eI))
158 rescale_velocities(ekind, cTC, 0, homenr, state->v);
163 // Set the T scaling lambda to 1 to have no scaling
164 // TODO: Do we have to do it on every non-t-couple step?
165 for (int i = 0; (i < inputrec->opts.ngtc); i++)
167 ekind->tcstat[i].lambda = 1.0;
172 void update_pcouple_before_coordinates(FILE* fplog,
174 const t_inputrec* inputrec,
176 matrix parrinellorahmanMu,
180 /* Berendsen P-coupling is completely handled after the coordinate update.
181 * Trotter P-coupling is handled by separate calls to trotter_update().
183 if (inputrec->epc == PressureCoupling::ParrinelloRahman
184 && do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
186 real dtpc = inputrec->nstpcouple * inputrec->delta_t;
188 parrinellorahman_pcoupl(
189 fplog, step, inputrec, dtpc, state->pres_prev, state->box, state->box_rel, state->boxv, M, parrinellorahmanMu, bInitStep);
193 void update_pcouple_after_coordinates(FILE* fplog,
195 const t_inputrec* inputrec,
197 gmx::ArrayRef<const unsigned short> cFREEZE,
198 const matrix pressure,
199 const matrix forceVirial,
200 const matrix constraintVirial,
201 matrix pressureCouplingMu,
204 gmx::BoxDeformation* boxDeformation,
205 const bool scaleCoordinates)
209 /* Cast to real for faster code, no loss in precision (see comment above) */
210 real dt = inputrec->delta_t;
213 /* now update boxes */
214 switch (inputrec->epc)
216 case (PressureCoupling::No): break;
217 case (PressureCoupling::Berendsen):
218 if (do_per_step(step, inputrec->nstpcouple))
220 real dtpc = inputrec->nstpcouple * dt;
221 pressureCouplingCalculateScalingMatrix<PressureCoupling::Berendsen>(fplog,
230 &state->baros_integral);
231 pressureCouplingScaleBoxAndCoordinates<PressureCoupling::Berendsen>(
239 gmx::ArrayRef<gmx::RVec>(),
245 case (PressureCoupling::CRescale):
246 if (do_per_step(step, inputrec->nstpcouple))
248 real dtpc = inputrec->nstpcouple * dt;
249 pressureCouplingCalculateScalingMatrix<PressureCoupling::CRescale>(fplog,
258 &state->baros_integral);
259 pressureCouplingScaleBoxAndCoordinates<PressureCoupling::CRescale>(inputrec,
272 case (PressureCoupling::ParrinelloRahman):
273 if (do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
275 /* The box velocities were updated in do_pr_pcoupl,
276 * but we dont change the box vectors until we get here
277 * since we need to be able to shift/unshift above.
279 real dtpc = inputrec->nstpcouple * dt;
280 for (int i = 0; i < DIM; i++)
282 for (int m = 0; m <= i; m++)
284 state->box[i][m] += dtpc * state->boxv[i][m];
287 preserve_box_shape(inputrec, state->box_rel, state->box);
289 /* Scale the coordinates */
290 if (scaleCoordinates)
292 auto* x = state->x.rvec_array();
293 for (int n = start; n < start + homenr; n++)
295 tmvmul_ur0(pressureCouplingMu, x[n], x[n]);
300 case (PressureCoupling::Mttk):
301 switch (inputrec->epct)
303 case (PressureCouplingType::Isotropic):
304 /* DIM * eta = ln V. so DIM*eta_new = DIM*eta_old + DIM*dt*veta =>
305 ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) =>
306 Side length scales as exp(veta*dt) */
308 msmul(state->box, std::exp(state->veta * dt), state->box);
310 /* Relate veta to boxv. veta = d(eta)/dT = (1/DIM)*1/V dV/dT.
311 o If we assume isotropic scaling, and box length scaling
312 factor L, then V = L^DIM (det(M)). So dV/dt = DIM
313 L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt. The
314 determinant of B is L^DIM det(M), and the determinant
315 of dB/dt is (dL/dT)^DIM det (M). veta will be
316 (det(dB/dT)/det(B))^(1/3). Then since M =
317 B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */
319 msmul(state->box, state->veta, state->boxv);
329 auto localX = makeArrayRef(state->x).subArray(start, homenr);
330 boxDeformation->apply(localX, state->box, step);
334 extern bool update_randomize_velocities(const t_inputrec* ir,
338 gmx::ArrayRef<const unsigned short> cTC,
339 gmx::ArrayRef<const real> invMass,
340 gmx::ArrayRef<gmx::RVec> v,
341 const gmx::Update* upd,
342 const gmx::Constraints* constr)
345 real rate = (ir->delta_t) / ir->opts.tau_t[0];
347 if (ir->etc == TemperatureCoupling::Andersen && constr != nullptr)
349 /* Currently, Andersen thermostat does not support constrained
350 systems. Functionality exists in the andersen_tcoupl
351 function in GROMACS 4.5.7 to allow this combination. That
352 code could be ported to the current random-number
353 generation approach, but has not yet been done because of
354 lack of time and resources. */
356 "Normal Andersen is currently not supported with constraints, use massive "
360 /* proceed with andersen if 1) it's fixed probability per
361 particle andersen or 2) it's massive andersen and it's tau_t/dt */
362 if ((ir->etc == TemperatureCoupling::Andersen) || do_per_step(step, gmx::roundToInt(1.0 / rate)))
365 ir, step, cr, homenr, cTC, invMass, v, rate, upd->getAndersenRandomizeGroup(), upd->getBoltzmanFactor());
371 /* these integration routines are only referenced inside this file */
372 static void NHC_trotter(const t_grpopts* opts,
374 const gmx_ekindata_t* ekind,
380 const t_extmass* MassQ,
384 /* general routine for both barostat and thermostat nose hoover chains */
387 double Ekin, Efac, reft, kT, nd;
392 int mstepsi, mstepsj;
393 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
394 int nh = opts->nhchainlength;
397 mstepsi = mstepsj = ns;
399 /* if scalefac is NULL, we are doing the NHC of the barostat */
402 if (scalefac == nullptr)
407 for (i = 0; i < nvar; i++)
410 /* make it easier to iterate by selecting
411 out the sub-array that corresponds to this T group */
415 gmx::ArrayRef<const double> iQinv;
418 iQinv = gmx::arrayRefFromArray(&MassQ->QPinv[i * nh], nh);
419 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
420 reft = std::max<real>(0, opts->ref_t[0]);
421 Ekin = gmx::square(*veta) / MassQ->Winv;
425 iQinv = gmx::arrayRefFromArray(&MassQ->Qinv[i * nh], nh);
426 const t_grp_tcstat* tcstat = &ekind->tcstat[i];
428 reft = std::max<real>(0, opts->ref_t[i]);
431 Ekin = 2 * trace(tcstat->ekinf) * tcstat->ekinscalef_nhc;
435 Ekin = 2 * trace(tcstat->ekinh) * tcstat->ekinscaleh_nhc;
438 kT = gmx::c_boltz * reft;
440 for (mi = 0; mi < mstepsi; mi++)
442 for (mj = 0; mj < mstepsj; mj++)
444 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
445 dt = sy_const[ns][mj] * dtfull / mstepsi;
447 /* compute the thermal forces */
448 GQ[0] = iQinv[0] * (Ekin - nd * kT);
450 for (j = 0; j < nh - 1; j++)
452 if (iQinv[j + 1] > 0)
454 /* we actually don't need to update here if we save the
455 state of the GQ, but it's easier to just recompute*/
456 GQ[j + 1] = iQinv[j + 1] * ((gmx::square(ivxi[j]) / iQinv[j]) - kT);
464 ivxi[nh - 1] += 0.25 * dt * GQ[nh - 1];
465 for (j = nh - 1; j > 0; j--)
467 Efac = exp(-0.125 * dt * ivxi[j]);
468 ivxi[j - 1] = Efac * (ivxi[j - 1] * Efac + 0.25 * dt * GQ[j - 1]);
471 Efac = exp(-0.5 * dt * ivxi[0]);
480 Ekin *= (Efac * Efac);
482 /* Issue - if the KE is an average of the last and the current temperatures, then we
483 might not be able to scale the kinetic energy directly with this factor. Might
484 take more bookkeeping -- have to think about this a bit more . . . */
486 GQ[0] = iQinv[0] * (Ekin - nd * kT);
488 /* update thermostat positions */
489 for (j = 0; j < nh; j++)
491 ixi[j] += 0.5 * dt * ivxi[j];
494 for (j = 0; j < nh - 1; j++)
496 Efac = exp(-0.125 * dt * ivxi[j + 1]);
497 ivxi[j] = Efac * (ivxi[j] * Efac + 0.25 * dt * GQ[j]);
498 if (iQinv[j + 1] > 0)
500 GQ[j + 1] = iQinv[j + 1] * ((gmx::square(ivxi[j]) / iQinv[j]) - kT);
507 ivxi[nh - 1] += 0.25 * dt * GQ[nh - 1];
514 static void boxv_trotter(const t_inputrec* ir,
518 const gmx_ekindata_t* ekind,
521 const t_extmass* MassQ)
528 tensor ekinmod, localpres;
530 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
531 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
534 if (ir->epct == PressureCouplingType::SemiIsotropic)
543 /* eta is in pure units. veta is in units of ps^-1. GW is in
544 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
545 taken to use only RATIOS of eta in updating the volume. */
547 /* we take the partial pressure tensors, modify the
548 kinetic energy tensor, and recovert to pressure */
550 if (ir->opts.nrdf[0] == 0)
552 gmx_fatal(FARGS, "Barostat is coupled to a T-group with no degrees of freedom\n");
554 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
555 alpha = 1.0 + DIM / (static_cast<double>(ir->opts.nrdf[0]));
556 alpha *= ekind->tcstat[0].ekinscalef_nhc;
557 msmul(ekind->ekin, alpha, ekinmod);
558 /* for now, we use Elr = 0, because if you want to get it right, you
559 really should be using PME. Maybe print a warning? */
561 pscal = calc_pres(ir->pbcType, nwall, box, ekinmod, vir, localpres) + pcorr;
564 GW = (vol * (MassQ->Winv / gmx::c_presfac))
565 * (DIM * pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
567 *veta += 0.5 * dt * GW;
571 * This file implements temperature and pressure coupling algorithms:
572 * For now only the Weak coupling and the modified weak coupling.
574 * Furthermore computation of pressure and temperature is done here
578 real calc_pres(PbcType pbcType, int nwall, const matrix box, const tensor ekin, const tensor vir, tensor pres)
583 if (pbcType == PbcType::No || (pbcType == PbcType::XY && nwall != 2))
589 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
590 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
594 fac = gmx::c_presfac * 2.0 / det(box);
595 for (n = 0; (n < DIM); n++)
597 for (m = 0; (m < DIM); m++)
599 pres[n][m] = (ekin[n][m] - vir[n][m]) * fac;
605 pr_rvecs(debug, 0, "PC: pres", pres, DIM);
606 pr_rvecs(debug, 0, "PC: ekin", ekin, DIM);
607 pr_rvecs(debug, 0, "PC: vir ", vir, DIM);
608 pr_rvecs(debug, 0, "PC: box ", box, DIM);
611 return trace(pres) / DIM;
614 real calc_temp(real ekin, real nrdf)
618 return (2.0 * ekin) / (nrdf * gmx::c_boltz);
626 /*! \brief Sets 1/mass for Parrinello-Rahman in wInv; NOTE: c_presfac is not included, so not in GROMACS units! */
627 static void calcParrinelloRahmanInvMass(const t_inputrec* ir, const matrix box, tensor wInv)
631 /* TODO: See if we can make the mass independent of the box size */
632 maxBoxLength = std::max(box[XX][XX], box[YY][YY]);
633 maxBoxLength = std::max(maxBoxLength, box[ZZ][ZZ]);
635 for (int d = 0; d < DIM; d++)
637 for (int n = 0; n < DIM; n++)
639 wInv[d][n] = (4 * M_PI * M_PI * ir->compress[d][n]) / (3 * ir->tau_p * ir->tau_p * maxBoxLength);
644 void parrinellorahman_pcoupl(FILE* fplog,
646 const t_inputrec* ir,
656 /* This doesn't do any coordinate updating. It just
657 * integrates the box vector equations from the calculated
658 * acceleration due to pressure difference. We also compute
659 * the tensor M which is used in update to couple the particle
660 * coordinates to the box vectors.
662 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
665 * M_nk = (h') * (h' * h + h' h) * h
667 * with the dots denoting time derivatives and h is the transformation from
668 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
669 * This also goes for the pressure and M tensors - they are transposed relative
670 * to ours. Our equation thus becomes:
673 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
675 * where b is the gromacs box matrix.
676 * Our box accelerations are given by
678 * b = vol/W inv(box') * (P-ref_P) (=h')
681 real vol = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
682 real atot, arel, change;
683 tensor invbox, pdiff, t1, t2;
685 gmx::invertBoxMatrix(box, invbox);
689 /* Note that c_presfac does not occur here.
690 * The pressure and compressibility always occur as a product,
691 * therefore the pressure unit drops out.
694 calcParrinelloRahmanInvMass(ir, box, winv);
696 m_sub(pres, ir->ref_p, pdiff);
698 if (ir->epct == PressureCouplingType::SurfaceTension)
700 /* Unlike Berendsen coupling it might not be trivial to include a z
701 * pressure correction here? On the other hand we don't scale the
702 * box momentarily, but change accelerations, so it might not be crucial.
704 real xy_pressure = 0.5 * (pres[XX][XX] + pres[YY][YY]);
705 for (int d = 0; d < ZZ; d++)
707 pdiff[d][d] = (xy_pressure - (pres[ZZ][ZZ] - ir->ref_p[d][d] / box[d][d]));
711 tmmul(invbox, pdiff, t1);
712 /* Move the off-diagonal elements of the 'force' to one side to ensure
713 * that we obey the box constraints.
715 for (int d = 0; d < DIM; d++)
717 for (int n = 0; n < d; n++)
719 t1[d][n] += t1[n][d];
726 case PressureCouplingType::Anisotropic:
727 for (int d = 0; d < DIM; d++)
729 for (int n = 0; n <= d; n++)
731 t1[d][n] *= winv[d][n] * vol;
735 case PressureCouplingType::Isotropic:
736 /* calculate total volume acceleration */
737 atot = box[XX][XX] * box[YY][YY] * t1[ZZ][ZZ] + box[XX][XX] * t1[YY][YY] * box[ZZ][ZZ]
738 + t1[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
739 arel = atot / (3 * vol);
740 /* set all RELATIVE box accelerations equal, and maintain total V
742 for (int d = 0; d < DIM; d++)
744 for (int n = 0; n <= d; n++)
746 t1[d][n] = winv[0][0] * vol * arel * box[d][n];
750 case PressureCouplingType::SemiIsotropic:
751 case PressureCouplingType::SurfaceTension:
752 /* Note the correction to pdiff above for surftens. coupling */
754 /* calculate total XY volume acceleration */
755 atot = box[XX][XX] * t1[YY][YY] + t1[XX][XX] * box[YY][YY];
756 arel = atot / (2 * box[XX][XX] * box[YY][YY]);
757 /* set RELATIVE XY box accelerations equal, and maintain total V
758 * change speed. Dont change the third box vector accelerations */
759 for (int d = 0; d < ZZ; d++)
761 for (int n = 0; n <= d; n++)
763 t1[d][n] = winv[d][n] * vol * arel * box[d][n];
766 for (int n = 0; n < DIM; n++)
768 t1[ZZ][n] *= winv[ZZ][n] * vol;
773 "Parrinello-Rahman pressure coupling type %s "
774 "not supported yet\n",
775 enumValueToString(ir->epct));
779 for (int d = 0; d < DIM; d++)
781 for (int n = 0; n <= d; n++)
783 boxv[d][n] += dt * t1[d][n];
785 /* We do NOT update the box vectors themselves here, since
786 * we need them for shifting later. It is instead done last
787 * in the update() routine.
790 /* Calculate the change relative to diagonal elements-
791 since it's perfectly ok for the off-diagonal ones to
792 be zero it doesn't make sense to check the change relative
796 change = std::fabs(dt * boxv[d][n] / box[d][d]);
798 if (change > maxchange)
805 if (maxchange > 0.01 && fplog)
809 "\nStep %s Warning: Pressure scaling more than 1%%. "
810 "This may mean your system\n is not yet equilibrated. "
811 "Use of Parrinello-Rahman pressure coupling during\n"
812 "equilibration can lead to simulation instability, "
813 "and is discouraged.\n",
814 gmx_step_str(step, buf));
818 preserve_box_shape(ir, box_rel, boxv);
820 mtmul(boxv, box, t1); /* t1=boxv * b' */
821 mmul(invbox, t1, t2);
822 mtmul(t2, invbox, M);
824 /* Determine the scaling matrix mu for the coordinates */
825 for (int d = 0; d < DIM; d++)
827 for (int n = 0; n <= d; n++)
829 t1[d][n] = box[d][n] + dt * boxv[d][n];
832 preserve_box_shape(ir, box_rel, t1);
833 /* t1 is the box at t+dt, determine mu as the relative change */
834 mmul_ur0(invbox, t1, mu);
837 //! Return compressibility factor for entry (i,j) of Berendsen / C-rescale scaling matrix
838 static inline real compressibilityFactor(int i, int j, const t_inputrec* ir, real dt)
840 return ir->compress[i][j] * dt / ir->tau_p;
843 //! Details of Berendsen / C-rescale scaling matrix calculation
844 template<PressureCoupling pressureCouplingType>
845 static void calculateScalingMatrixImplDetail(const t_inputrec* ir,
850 real scalar_pressure,
854 //! Calculate Berendsen / C-rescale scaling matrix
855 template<PressureCoupling pressureCouplingType>
856 static void calculateScalingMatrixImpl(const t_inputrec* ir,
863 real scalar_pressure = 0;
864 real xy_pressure = 0;
865 for (int d = 0; d < DIM; d++)
867 scalar_pressure += pres[d][d] / DIM;
870 xy_pressure += pres[d][d] / (DIM - 1);
874 calculateScalingMatrixImplDetail<pressureCouplingType>(
875 ir, mu, dt, pres, box, scalar_pressure, xy_pressure, step);
879 void calculateScalingMatrixImplDetail<PressureCoupling::Berendsen>(const t_inputrec* ir,
884 real scalar_pressure,
886 int64_t gmx_unused step)
891 case PressureCouplingType::Isotropic:
892 for (int d = 0; d < DIM; d++)
894 mu[d][d] = 1.0 - compressibilityFactor(d, d, ir, dt) * (ir->ref_p[d][d] - scalar_pressure) / DIM;
897 case PressureCouplingType::SemiIsotropic:
898 for (int d = 0; d < ZZ; d++)
900 mu[d][d] = 1.0 - compressibilityFactor(d, d, ir, dt) * (ir->ref_p[d][d] - xy_pressure) / DIM;
903 1.0 - compressibilityFactor(ZZ, ZZ, ir, dt) * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM;
905 case PressureCouplingType::Anisotropic:
906 for (int d = 0; d < DIM; d++)
908 for (int n = 0; n < DIM; n++)
910 mu[d][n] = (d == n ? 1.0 : 0.0)
911 - compressibilityFactor(d, n, ir, dt) * (ir->ref_p[d][n] - pres[d][n]) / DIM;
915 case PressureCouplingType::SurfaceTension:
916 /* ir->ref_p[0/1] is the reference surface-tension times *
917 * the number of surfaces */
918 if (ir->compress[ZZ][ZZ] != 0.0F)
920 p_corr_z = dt / ir->tau_p * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
924 /* when the compressibity is zero, set the pressure correction *
925 * in the z-direction to zero to get the correct surface tension */
928 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ] * p_corr_z;
929 for (int d = 0; d < DIM - 1; d++)
932 + compressibilityFactor(d, d, ir, dt)
933 * (ir->ref_p[d][d] / (mu[ZZ][ZZ] * box[ZZ][ZZ])
934 - (pres[ZZ][ZZ] + p_corr_z - xy_pressure))
940 "Berendsen pressure coupling type %s not supported yet\n",
941 enumValueToString(ir->epct));
946 void calculateScalingMatrixImplDetail<PressureCoupling::CRescale>(const t_inputrec* ir,
951 real scalar_pressure,
955 gmx::ThreeFry2x64<64> rng(ir->ld_seed, gmx::RandomDomain::Barostat);
956 gmx::NormalDistribution<real> normalDist;
957 rng.restart(step, 0);
959 for (int d = 0; d < DIM; d++)
965 real kt = ir->opts.ref_t[0] * gmx::c_boltz;
973 case PressureCouplingType::Isotropic:
974 gauss = normalDist(rng);
975 for (int d = 0; d < DIM; d++)
977 const real factor = compressibilityFactor(d, d, ir, dt);
978 mu[d][d] = std::exp(-factor * (ir->ref_p[d][d] - scalar_pressure) / DIM
979 + std::sqrt(2.0 * kt * factor * gmx::c_presfac / vol) * gauss / DIM);
982 case PressureCouplingType::SemiIsotropic:
983 gauss = normalDist(rng);
984 gauss2 = normalDist(rng);
985 for (int d = 0; d < ZZ; d++)
987 const real factor = compressibilityFactor(d, d, ir, dt);
988 mu[d][d] = std::exp(-factor * (ir->ref_p[d][d] - xy_pressure) / DIM
989 + std::sqrt((DIM - 1) * 2.0 * kt * factor * gmx::c_presfac / vol / DIM)
990 / (DIM - 1) * gauss);
993 const real factor = compressibilityFactor(ZZ, ZZ, ir, dt);
994 mu[ZZ][ZZ] = std::exp(-factor * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM
995 + std::sqrt(2.0 * kt * factor * gmx::c_presfac / vol / DIM) * gauss2);
998 case PressureCouplingType::SurfaceTension:
999 gauss = normalDist(rng);
1000 gauss2 = normalDist(rng);
1001 for (int d = 0; d < ZZ; d++)
1003 const real factor = compressibilityFactor(d, d, ir, dt);
1004 /* Notice: we here use ref_p[ZZ][ZZ] as isotropic pressure and ir->ref_p[d][d] as surface tension */
1005 mu[d][d] = std::exp(
1006 -factor * (ir->ref_p[ZZ][ZZ] - ir->ref_p[d][d] / box[ZZ][ZZ] - xy_pressure) / DIM
1007 + std::sqrt(4.0 / 3.0 * kt * factor * gmx::c_presfac / vol) / (DIM - 1) * gauss);
1010 const real factor = compressibilityFactor(ZZ, ZZ, ir, dt);
1011 mu[ZZ][ZZ] = std::exp(-factor * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM
1012 + std::sqrt(2.0 / 3.0 * kt * factor * gmx::c_presfac / vol) * gauss2);
1017 "C-rescale pressure coupling type %s not supported yet\n",
1018 enumValueToString(ir->epct));
1022 template<PressureCoupling pressureCouplingType>
1023 void pressureCouplingCalculateScalingMatrix(FILE* fplog,
1025 const t_inputrec* ir,
1029 const matrix force_vir,
1030 const matrix constraint_vir,
1032 double* baros_integral)
1034 // If the support here increases, we need to also add the template declarations
1035 // for new cases below.
1036 static_assert(pressureCouplingType == PressureCoupling::Berendsen
1037 || pressureCouplingType == PressureCoupling::CRescale,
1038 "pressureCouplingCalculateScalingMatrix is only implemented for Berendsen and "
1039 "C-rescale pressure coupling");
1041 calculateScalingMatrixImpl<pressureCouplingType>(ir, mu, dt, pres, box, step);
1043 /* To fulfill the orientation restrictions on triclinic boxes
1044 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
1045 * the other elements of mu to first order.
1047 mu[YY][XX] += mu[XX][YY];
1048 mu[ZZ][XX] += mu[XX][ZZ];
1049 mu[ZZ][YY] += mu[YY][ZZ];
1054 /* Keep track of the work the barostat applies on the system.
1055 * Without constraints force_vir tells us how Epot changes when scaling.
1056 * With constraints constraint_vir gives us the constraint contribution
1057 * to both Epot and Ekin. Although we are not scaling velocities, scaling
1058 * the coordinates leads to scaling of distances involved in constraints.
1059 * This in turn changes the angular momentum (even if the constrained
1060 * distances are corrected at the next step). The kinetic component
1061 * of the constraint virial captures the angular momentum change.
1063 for (int d = 0; d < DIM; d++)
1065 for (int n = 0; n <= d; n++)
1068 2 * (mu[d][n] - (n == d ? 1 : 0)) * (force_vir[d][n] + constraint_vir[d][n]);
1074 pr_rvecs(debug, 0, "PC: pres ", pres, 3);
1075 pr_rvecs(debug, 0, "PC: mu ", mu, 3);
1078 if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 || mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01
1079 || mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
1084 "\nStep %s Warning: pressure scaling more than 1%%, "
1086 gmx_step_str(step, buf2),
1092 fprintf(fplog, "%s", buf);
1094 fprintf(stderr, "%s", buf);
1098 template void pressureCouplingCalculateScalingMatrix<PressureCoupling::CRescale>(FILE*,
1109 template void pressureCouplingCalculateScalingMatrix<PressureCoupling::Berendsen>(FILE*,
1120 template<PressureCoupling pressureCouplingType>
1121 void pressureCouplingScaleBoxAndCoordinates(const t_inputrec* ir,
1127 gmx::ArrayRef<gmx::RVec> x,
1128 gmx::ArrayRef<gmx::RVec> v,
1129 gmx::ArrayRef<const unsigned short> cFREEZE,
1131 const bool scaleCoordinates)
1133 // If the support here increases, we need to also add the template declarations
1134 // for new cases below.
1135 static_assert(pressureCouplingType == PressureCoupling::Berendsen
1136 || pressureCouplingType == PressureCoupling::CRescale,
1137 "pressureCouplingScaleBoxAndCoordinates is only implemented for Berendsen and "
1138 "C-rescale pressure coupling");
1140 ivec* nFreeze = ir->opts.nFreeze;
1142 if (pressureCouplingType == PressureCoupling::CRescale)
1144 gmx::invertBoxMatrix(mu, inv_mu);
1147 /* Scale the positions and the velocities */
1148 if (scaleCoordinates)
1150 const int gmx_unused numThreads = gmx_omp_nthreads_get(ModuleMultiThread::Update);
1151 #pragma omp parallel for num_threads(numThreads) schedule(static)
1152 for (int n = start; n < start + nr_atoms; n++)
1154 // Trivial OpenMP region that does not throw
1156 if (!cFREEZE.empty())
1161 if (!nFreeze[g][XX])
1163 x[n][XX] = mu[XX][XX] * x[n][XX] + mu[YY][XX] * x[n][YY] + mu[ZZ][XX] * x[n][ZZ];
1164 if (pressureCouplingType == PressureCoupling::CRescale)
1166 v[n][XX] = inv_mu[XX][XX] * v[n][XX] + inv_mu[YY][XX] * v[n][YY]
1167 + inv_mu[ZZ][XX] * v[n][ZZ];
1170 if (!nFreeze[g][YY])
1172 x[n][YY] = mu[YY][YY] * x[n][YY] + mu[ZZ][YY] * x[n][ZZ];
1173 if (pressureCouplingType == PressureCoupling::CRescale)
1175 v[n][YY] = inv_mu[YY][YY] * v[n][YY] + inv_mu[ZZ][YY] * v[n][ZZ];
1178 if (!nFreeze[g][ZZ])
1180 x[n][ZZ] = mu[ZZ][ZZ] * x[n][ZZ];
1181 if (pressureCouplingType == PressureCoupling::CRescale)
1183 v[n][ZZ] = inv_mu[ZZ][ZZ] * v[n][ZZ];
1188 /* compute final boxlengths */
1189 for (int d = 0; d < DIM; d++)
1191 box[d][XX] = mu[XX][XX] * box[d][XX] + mu[YY][XX] * box[d][YY] + mu[ZZ][XX] * box[d][ZZ];
1192 box[d][YY] = mu[YY][YY] * box[d][YY] + mu[ZZ][YY] * box[d][ZZ];
1193 box[d][ZZ] = mu[ZZ][ZZ] * box[d][ZZ];
1196 preserve_box_shape(ir, box_rel, box);
1198 /* (un)shifting should NOT be done after this,
1199 * since the box vectors might have changed
1201 inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
1205 pressureCouplingScaleBoxAndCoordinates<PressureCoupling::Berendsen>(const t_inputrec*,
1211 gmx::ArrayRef<gmx::RVec>,
1212 gmx::ArrayRef<gmx::RVec>,
1213 gmx::ArrayRef<const unsigned short>,
1219 pressureCouplingScaleBoxAndCoordinates<PressureCoupling::CRescale>(const t_inputrec*,
1225 gmx::ArrayRef<gmx::RVec>,
1226 gmx::ArrayRef<gmx::RVec>,
1227 gmx::ArrayRef<const unsigned short>,
1231 void berendsen_tcoupl(const t_inputrec* ir, gmx_ekindata_t* ekind, real dt, std::vector<double>& therm_integral)
1233 const t_grpopts* opts = &ir->opts;
1235 for (int i = 0; (i < opts->ngtc); i++)
1239 if (ir->eI == IntegrationAlgorithm::VV)
1241 Ek = trace(ekind->tcstat[i].ekinf);
1242 T = ekind->tcstat[i].T;
1246 Ek = trace(ekind->tcstat[i].ekinh);
1247 T = ekind->tcstat[i].Th;
1250 if ((opts->tau_t[i] > 0) && (T > 0.0))
1252 real reft = std::max<real>(0, opts->ref_t[i]);
1253 real lll = std::sqrt(1.0 + (dt / opts->tau_t[i]) * (reft / T - 1.0));
1254 ekind->tcstat[i].lambda = std::max<real>(std::min<real>(lll, 1.25), 0.8);
1258 ekind->tcstat[i].lambda = 1.0;
1261 /* Keep track of the amount of energy we are adding to the system */
1262 therm_integral[i] -= (gmx::square(ekind->tcstat[i].lambda) - 1) * Ek;
1266 fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n", i, T, ekind->tcstat[i].lambda);
1271 void andersen_tcoupl(const t_inputrec* ir,
1273 const t_commrec* cr,
1275 gmx::ArrayRef<const unsigned short> cTC,
1276 gmx::ArrayRef<const real> invMass,
1277 gmx::ArrayRef<gmx::RVec> v,
1279 const std::vector<bool>& randomize,
1280 gmx::ArrayRef<const real> boltzfac)
1282 const int* gatindex = (DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1285 gmx::ThreeFry2x64<0> rng(ir->andersen_seed, gmx::RandomDomain::Thermostat);
1286 gmx::UniformRealDistribution<real> uniformDist;
1287 gmx::TabulatedNormalDistribution<real, 14> normalDist;
1289 /* randomize the velocities of the selected particles */
1291 for (i = 0; i < homenr; i++) /* now loop over the list of atoms */
1293 int ng = gatindex ? gatindex[i] : i;
1296 rng.restart(step, ng);
1300 gc = cTC[i]; /* assign the atom to a temperature group if there are more than one */
1304 if (ir->etc == TemperatureCoupling::AndersenMassive)
1306 /* Randomize particle always */
1311 /* Randomize particle probabilistically */
1312 uniformDist.reset();
1313 bRandomize = uniformDist(rng) < rate;
1320 scal = std::sqrt(boltzfac[gc] * invMass[i]);
1324 for (d = 0; d < DIM; d++)
1326 v[i][d] = scal * normalDist(rng);
1334 void nosehoover_tcoupl(const t_grpopts* opts,
1335 const gmx_ekindata_t* ekind,
1337 gmx::ArrayRef<double> xi,
1338 gmx::ArrayRef<double> vxi,
1339 const t_extmass* MassQ)
1344 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
1346 for (i = 0; (i < opts->ngtc); i++)
1348 reft = std::max<real>(0, opts->ref_t[i]);
1350 vxi[i] += dt * MassQ->Qinv[i] * (ekind->tcstat[i].Th - reft);
1351 xi[i] += dt * (oldvxi + vxi[i]) * 0.5;
1355 void trotter_update(const t_inputrec* ir,
1357 gmx_ekindata_t* ekind,
1358 const gmx_enerdata_t* enerd,
1362 gmx::ArrayRef<const unsigned short> cTC,
1363 gmx::ArrayRef<const real> invMass,
1364 const t_extmass* MassQ,
1365 gmx::ArrayRef<std::vector<int>> trotter_seqlist,
1366 TrotterSequence trotter_seqno)
1369 int n, i, d, ngtc, gc = 0, t;
1370 t_grp_tcstat* tcstat;
1371 const t_grpopts* opts;
1374 double * scalefac, dtc;
1375 rvec sumv = { 0, 0, 0 };
1378 if (trotter_seqno <= TrotterSequence::Two)
1380 step_eff = step - 1; /* the velocity verlet calls are actually out of order -- the first
1381 half step is actually the last half step from the previous step.
1382 Thus the first half step actually corresponds to the n-1 step*/
1389 bCouple = (ir->nsttcouple == 1 || do_per_step(step_eff + ir->nsttcouple, ir->nsttcouple));
1391 const gmx::ArrayRef<const int> trotter_seq = trotter_seqlist[static_cast<int>(trotter_seqno)];
1393 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
1397 dtc = ir->nsttcouple * ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
1398 opts = &(ir->opts); /* just for ease of referencing */
1401 snew(scalefac, opts->ngtc);
1402 for (i = 0; i < ngtc; i++)
1406 /* execute the series of trotter updates specified in the trotterpart array */
1408 for (i = 0; i < NTROTTERPARTS; i++)
1410 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
1411 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2)
1412 || (trotter_seq[i] == etrtNHC2))
1421 auto v = makeArrayRef(state->v);
1422 switch (trotter_seq[i])
1426 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir, enerd->term[F_PDISPCORR], MassQ);
1434 state->nhpres_xi.data(),
1435 state->nhpres_vxi.data(),
1447 state->nosehoover_xi.data(),
1448 state->nosehoover_vxi.data(),
1452 (ir->eI == IntegrationAlgorithm::VV));
1453 /* need to rescale the kinetic energies and velocities here. Could
1454 scale the velocities later, but we need them scaled in order to
1455 produce the correct outputs, so we'll scale them here. */
1457 for (t = 0; t < ngtc; t++)
1459 tcstat = &ekind->tcstat[t];
1460 tcstat->vscale_nhc = scalefac[t];
1461 tcstat->ekinscaleh_nhc *= (scalefac[t] * scalefac[t]);
1462 tcstat->ekinscalef_nhc *= (scalefac[t] * scalefac[t]);
1464 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
1465 /* but do we actually need the total? */
1467 /* modify the velocities as well */
1468 for (n = 0; n < homenr; n++)
1470 if (!cTC.empty()) /* does this conditional need to be here? is this always true?*/
1474 for (d = 0; d < DIM; d++)
1476 v[n][d] *= scalefac[gc];
1481 for (d = 0; d < DIM; d++)
1483 sumv[d] += (v[n][d]) / invMass[n];
1491 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
1496 extern void init_npt_masses(const t_inputrec* ir, t_state* state, t_extmass* MassQ, bool bInit)
1498 int n, i, j, d, ngtc, nh;
1499 const t_grpopts* opts;
1500 real reft, kT, ndj, nd;
1502 opts = &(ir->opts); /* just for ease of referencing */
1503 ngtc = ir->opts.ngtc;
1504 nh = state->nhchainlength;
1506 if (ir->eI == IntegrationAlgorithm::MD)
1510 MassQ->Qinv.resize(ngtc);
1512 for (i = 0; (i < ngtc); i++)
1514 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1516 MassQ->Qinv[i] = 1.0 / (gmx::square(opts->tau_t[i] / M_2PI) * opts->ref_t[i]);
1520 MassQ->Qinv[i] = 0.0;
1524 else if (EI_VV(ir->eI))
1526 /* Set pressure variables */
1530 if (state->vol0 == 0)
1532 state->vol0 = det(state->box);
1533 /* because we start by defining a fixed
1534 compressibility, we need the volume at this
1535 compressibility to solve the problem. */
1539 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1540 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1541 MassQ->Winv = (gmx::c_presfac * trace(ir->compress) * gmx::c_boltz * opts->ref_t[0])
1542 / (DIM * state->vol0 * gmx::square(ir->tau_p / M_2PI));
1543 /* An alternate mass definition, from Tuckerman et al. */
1544 /* MassQ->Winv = 1.0/(gmx::square(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*c_boltz*opts->ref_t[0]); */
1545 for (d = 0; d < DIM; d++)
1547 for (n = 0; n < DIM; n++)
1549 MassQ->Winvm[d][n] = gmx::c_presfac * ir->compress[d][n]
1550 / (state->vol0 * gmx::square(ir->tau_p / M_2PI));
1551 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1552 before using MTTK for anisotropic states.*/
1555 /* Allocate space for thermostat variables */
1558 MassQ->Qinv.resize(ngtc * nh);
1561 /* now, set temperature variables */
1562 for (i = 0; i < ngtc; i++)
1564 if (opts->tau_t[i] > 0 && opts->ref_t[i] > 0 && opts->nrdf[i] > 0)
1566 reft = std::max<real>(0, opts->ref_t[i]);
1568 kT = gmx::c_boltz * reft;
1569 for (j = 0; j < nh; j++)
1579 MassQ->Qinv[i * nh + j] = 1.0 / (gmx::square(opts->tau_t[i] / M_2PI) * ndj * kT);
1584 for (j = 0; j < nh; j++)
1586 MassQ->Qinv[i * nh + j] = 0.0;
1593 gmx::EnumerationArray<TrotterSequence, std::vector<int>>
1594 init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* MassQ, bool bTrotter)
1596 int i, j, nnhpres, nh;
1597 const t_grpopts* opts;
1598 real bmass, qmass, reft, kT;
1600 opts = &(ir->opts); /* just for ease of referencing */
1601 nnhpres = state->nnhpres;
1602 nh = state->nhchainlength;
1604 if (EI_VV(ir->eI) && (ir->epc == PressureCoupling::Mttk) && (ir->etc != TemperatureCoupling::NoseHoover))
1606 gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
1609 init_npt_masses(ir, state, MassQ, TRUE);
1611 /* first, initialize clear all the trotter calls */
1612 gmx::EnumerationArray<TrotterSequence, std::vector<int>> trotter_seq;
1613 for (i = 0; i < static_cast<int>(TrotterSequence::Count); i++)
1615 trotter_seq[i].resize(NTROTTERPARTS, etrtNONE);
1616 trotter_seq[i][0] = etrtSKIPALL;
1621 /* no trotter calls, so we never use the values in the array.
1622 * We access them (so we need to define them, but ignore
1628 /* compute the kinetic energy by using the half step velocities or
1629 * the kinetic energies, depending on the order of the trotter calls */
1631 if (ir->eI == IntegrationAlgorithm::VV)
1633 if (inputrecNptTrotter(ir))
1635 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1636 We start with the initial one. */
1637 /* first, a round that estimates veta. */
1638 trotter_seq[0][0] = etrtBAROV;
1640 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1642 /* The first half trotter update */
1643 trotter_seq[2][0] = etrtBAROV;
1644 trotter_seq[2][1] = etrtNHC;
1645 trotter_seq[2][2] = etrtBARONHC;
1647 /* The second half trotter update */
1648 trotter_seq[3][0] = etrtBARONHC;
1649 trotter_seq[3][1] = etrtNHC;
1650 trotter_seq[3][2] = etrtBAROV;
1652 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1654 else if (inputrecNvtTrotter(ir))
1656 /* This is the easy version - there are only two calls, both the same.
1657 Otherwise, even easier -- no calls */
1658 trotter_seq[2][0] = etrtNHC;
1659 trotter_seq[3][0] = etrtNHC;
1661 else if (inputrecNphTrotter(ir))
1663 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1664 We start with the initial one. */
1665 /* first, a round that estimates veta. */
1666 trotter_seq[0][0] = etrtBAROV;
1668 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1670 /* The first half trotter update */
1671 trotter_seq[2][0] = etrtBAROV;
1672 trotter_seq[2][1] = etrtBARONHC;
1674 /* The second half trotter update */
1675 trotter_seq[3][0] = etrtBARONHC;
1676 trotter_seq[3][1] = etrtBAROV;
1678 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1681 else if (ir->eI == IntegrationAlgorithm::VVAK)
1683 if (inputrecNptTrotter(ir))
1685 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1686 We start with the initial one. */
1687 /* first, a round that estimates veta. */
1688 trotter_seq[0][0] = etrtBAROV;
1690 /* The first half trotter update, part 1 -- double update, because it commutes */
1691 trotter_seq[1][0] = etrtNHC;
1693 /* The first half trotter update, part 2 */
1694 trotter_seq[2][0] = etrtBAROV;
1695 trotter_seq[2][1] = etrtBARONHC;
1697 /* The second half trotter update, part 1 */
1698 trotter_seq[3][0] = etrtBARONHC;
1699 trotter_seq[3][1] = etrtBAROV;
1701 /* The second half trotter update */
1702 trotter_seq[4][0] = etrtNHC;
1704 else if (inputrecNvtTrotter(ir))
1706 /* This is the easy version - there is only one call, both the same.
1707 Otherwise, even easier -- no calls */
1708 trotter_seq[1][0] = etrtNHC;
1709 trotter_seq[4][0] = etrtNHC;
1711 else if (inputrecNphTrotter(ir))
1713 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1714 We start with the initial one. */
1715 /* first, a round that estimates veta. */
1716 trotter_seq[0][0] = etrtBAROV;
1718 /* The first half trotter update, part 1 -- leave zero */
1719 trotter_seq[1][0] = etrtNHC;
1721 /* The first half trotter update, part 2 */
1722 trotter_seq[2][0] = etrtBAROV;
1723 trotter_seq[2][1] = etrtBARONHC;
1725 /* The second half trotter update, part 1 */
1726 trotter_seq[3][0] = etrtBARONHC;
1727 trotter_seq[3][1] = etrtBAROV;
1729 /* The second half trotter update -- blank for now */
1735 case PressureCouplingType::Isotropic:
1736 default: bmass = DIM * DIM; /* recommended mass parameters for isotropic barostat */
1739 MassQ->QPinv.resize(nnhpres * opts->nhchainlength);
1741 /* barostat temperature */
1742 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1744 reft = std::max<real>(0, opts->ref_t[0]);
1745 kT = gmx::c_boltz * reft;
1746 for (i = 0; i < nnhpres; i++)
1748 for (j = 0; j < nh; j++)
1758 MassQ->QPinv[i * opts->nhchainlength + j] =
1759 1.0 / (gmx::square(opts->tau_t[0] / M_2PI) * qmass * kT);
1765 for (i = 0; i < nnhpres; i++)
1767 for (j = 0; j < nh; j++)
1769 MassQ->QPinv[i * nh + j] = 0.0;
1776 static real energyNoseHoover(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1780 int nh = state->nhchainlength;
1782 for (int i = 0; i < ir->opts.ngtc; i++)
1784 const double* ixi = &state->nosehoover_xi[i * nh];
1785 const double* ivxi = &state->nosehoover_vxi[i * nh];
1786 const double* iQinv = &(MassQ->Qinv[i * nh]);
1788 real nd = ir->opts.nrdf[i];
1789 real reft = std::max<real>(ir->opts.ref_t[i], 0);
1790 real kT = gmx::c_boltz * reft;
1794 if (inputrecNvtTrotter(ir) || inputrecNptTrotter(ir))
1796 /* contribution from the thermal momenta of the NH chain */
1797 for (int j = 0; j < nh; j++)
1801 energy += 0.5 * gmx::square(ivxi[j]) / iQinv[j];
1802 /* contribution from the thermal variable of the NH chain */
1812 energy += ndj * ixi[j] * kT;
1816 else /* Other non Trotter temperature NH control -- no chains yet. */
1818 energy += 0.5 * gmx::c_boltz * nd * gmx::square(ivxi[0]) / iQinv[0];
1819 energy += nd * ixi[0] * kT;
1827 /* Returns the energy from the barostat thermostat chain */
1828 static real energyPressureMTTK(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1832 int nh = state->nhchainlength;
1834 for (int i = 0; i < state->nnhpres; i++)
1836 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1837 real reft = std::max<real>(ir->opts.ref_t[0], 0.0); /* using 'System' temperature */
1838 real kT = gmx::c_boltz * reft;
1840 for (int j = 0; j < nh; j++)
1842 double iQinv = MassQ->QPinv[i * nh + j];
1845 energy += 0.5 * gmx::square(state->nhpres_vxi[i * nh + j]) / iQinv;
1846 /* contribution from the thermal variable of the NH chain */
1847 energy += state->nhpres_xi[i * nh + j] * kT;
1852 "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f",
1855 state->nhpres_vxi[i * nh + j],
1856 state->nhpres_xi[i * nh + j]);
1864 /* Returns the energy accumulated by the V-rescale or Berendsen thermostat */
1865 static real energyVrescale(const t_inputrec* ir, const t_state* state)
1868 for (int i = 0; i < ir->opts.ngtc; i++)
1870 energy += state->therm_integral[i];
1876 real NPT_energy(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1880 if (ir->epc != PressureCoupling::No)
1882 /* Compute the contribution of the pressure to the conserved quantity*/
1884 real vol = det(state->box);
1888 case PressureCoupling::ParrinelloRahman:
1890 /* contribution from the pressure momenta */
1892 calcParrinelloRahmanInvMass(ir, state->box, invMass);
1893 for (int d = 0; d < DIM; d++)
1895 for (int n = 0; n <= d; n++)
1897 if (invMass[d][n] > 0)
1899 energyNPT += 0.5 * gmx::square(state->boxv[d][n])
1900 / (invMass[d][n] * gmx::c_presfac);
1905 /* Contribution from the PV term.
1906 * Not that with non-zero off-diagonal reference pressures,
1907 * i.e. applied shear stresses, there are additional terms.
1908 * We don't support this here, since that requires keeping
1909 * track of unwrapped box diagonal elements. This case is
1910 * excluded in integratorHasConservedEnergyQuantity().
1912 energyNPT += vol * trace(ir->ref_p) / (DIM * gmx::c_presfac);
1915 case PressureCoupling::Mttk:
1916 /* contribution from the pressure momenta */
1917 energyNPT += 0.5 * gmx::square(state->veta) / MassQ->Winv;
1919 /* contribution from the PV term */
1920 energyNPT += vol * trace(ir->ref_p) / (DIM * gmx::c_presfac);
1922 if (ir->epc == PressureCoupling::Mttk)
1924 /* contribution from the MTTK chain */
1925 energyNPT += energyPressureMTTK(ir, state, MassQ);
1928 case PressureCoupling::Berendsen:
1929 case PressureCoupling::CRescale: energyNPT += state->baros_integral; break;
1933 "Conserved energy quantity for pressure coupling is not handled. A case "
1934 "should be added with either the conserved quantity added or nothing added "
1935 "and an exclusion added to integratorHasConservedEnergyQuantity().");
1941 case TemperatureCoupling::No: break;
1942 case TemperatureCoupling::VRescale:
1943 case TemperatureCoupling::Berendsen: energyNPT += energyVrescale(ir, state); break;
1944 case TemperatureCoupling::NoseHoover:
1945 energyNPT += energyNoseHoover(ir, state, MassQ);
1947 case TemperatureCoupling::Andersen:
1948 case TemperatureCoupling::AndersenMassive:
1949 // Not supported, excluded in integratorHasConservedEnergyQuantity()
1954 "Conserved energy quantity for temperature coupling is not handled. A case "
1955 "should be added with either the conserved quantity added or nothing added and "
1956 "an exclusion added to integratorHasConservedEnergyQuantity().");
1963 static real vrescale_sumnoises(real nn, gmx::ThreeFry2x64<>* rng, gmx::NormalDistribution<real>* normalDist)
1966 * Returns the sum of nn independent gaussian noises squared
1967 * (i.e. equivalent to summing the square of the return values
1968 * of nn calls to a normal distribution).
1970 const real ndeg_tol = 0.0001;
1972 gmx::GammaDistribution<real> gammaDist(0.5 * nn, 1.0);
1974 if (nn < 2 + ndeg_tol)
1979 nn_int = gmx::roundToInt(nn);
1981 if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1984 "The v-rescale thermostat was called with a group with #DOF=%f, but for "
1985 "#DOF<3 only integer #DOF are supported",
1990 for (i = 0; i < nn_int; i++)
1992 gauss = (*normalDist)(*rng);
1998 /* Use a gamma distribution for any real nn > 2 */
1999 r = 2.0 * gammaDist(*rng);
2005 real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut, int64_t step, int64_t seed)
2008 * Generates a new value for the kinetic energy,
2009 * according to Bussi et al JCP (2007), Eq. (A7)
2010 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
2011 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
2012 * ndeg: number of degrees of freedom of the atoms to be thermalized
2013 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
2015 real factor, rr, ekin_new;
2016 gmx::ThreeFry2x64<64> rng(seed, gmx::RandomDomain::Thermostat);
2017 gmx::NormalDistribution<real> normalDist;
2021 factor = exp(-1.0 / taut);
2028 rng.restart(step, 0);
2030 rr = normalDist(rng);
2034 * (sigma * (vrescale_sumnoises(ndeg - 1, &rng, &normalDist) + rr * rr) / ndeg - kk)
2035 + 2.0 * rr * std::sqrt(kk * sigma / ndeg * (1.0 - factor) * factor);
2040 void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, real dt, gmx::ArrayRef<double> therm_integral)
2042 const t_grpopts* opts;
2044 real Ek, Ek_ref1, Ek_ref, Ek_new;
2048 for (i = 0; (i < opts->ngtc); i++)
2050 if (ir->eI == IntegrationAlgorithm::VV)
2052 Ek = trace(ekind->tcstat[i].ekinf);
2056 Ek = trace(ekind->tcstat[i].ekinh);
2059 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
2061 Ek_ref1 = 0.5 * opts->ref_t[i] * gmx::c_boltz;
2062 Ek_ref = Ek_ref1 * opts->nrdf[i];
2064 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i], opts->tau_t[i] / dt, step, ir->ld_seed);
2066 /* Analytically Ek_new>=0, but we check for rounding errors */
2069 ekind->tcstat[i].lambda = 0.0;
2073 ekind->tcstat[i].lambda = std::sqrt(Ek_new / Ek);
2076 therm_integral[i] -= Ek_new - Ek;
2081 "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
2086 ekind->tcstat[i].lambda);
2091 ekind->tcstat[i].lambda = 1.0;
2096 void rescale_velocities(const gmx_ekindata_t* ekind,
2097 gmx::ArrayRef<const unsigned short> cTC,
2100 gmx::ArrayRef<gmx::RVec> v)
2102 gmx::ArrayRef<const t_grp_tcstat> tcstat = ekind->tcstat;
2104 for (int n = start; n < end; n++)
2111 const real lg = tcstat[gt].lambda;
2112 for (int d = 0; d < DIM; d++)
2119 //! Check whether we're doing simulated annealing
2120 bool doSimulatedAnnealing(const t_inputrec* ir)
2122 for (int i = 0; i < ir->opts.ngtc; i++)
2124 /* set bSimAnn if any group is being annealed */
2125 if (ir->opts.annealing[i] != SimulatedAnnealing::No)
2133 // TODO If we keep simulated annealing, make a proper module that
2134 // does not rely on changing inputrec.
2135 bool initSimulatedAnnealing(t_inputrec* ir, gmx::Update* upd)
2137 bool doSimAnnealing = doSimulatedAnnealing(ir);
2140 update_annealing_target_temp(ir, ir->init_t, upd);
2142 return doSimAnnealing;
2145 /* set target temperatures if we are annealing */
2146 void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd)
2148 for (int temperatureGroup = 0; temperatureGroup < ir->opts.ngtc; temperatureGroup++)
2150 const auto& inputrec = *ir;
2151 if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::No)
2156 inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Single
2157 || inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Periodic,
2158 gmx::formatString("Unknown simulated annealing algorithm for temperature group %d", temperatureGroup)
2161 const auto npoints = inputrec.opts.anneal_npoints[temperatureGroup];
2162 if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Periodic)
2164 /* calculate time modulo the period */
2165 const auto pert = inputrec.opts.anneal_time[temperatureGroup][npoints - 1];
2166 const auto n = static_cast<int>(t / pert);
2167 thist = t - n * pert; /* modulo time */
2168 /* Make sure rounding didn't get us outside the interval */
2169 if (std::fabs(thist - pert) < GMX_REAL_EPS * 100)
2174 else if (inputrec.opts.annealing[temperatureGroup] == SimulatedAnnealing::Single)
2178 /* We are doing annealing for this group if we got here,
2179 * and we have the (relative) time as thist.
2180 * calculate target temp */
2182 while ((j < npoints - 1) && (thist > (inputrec.opts.anneal_time[temperatureGroup][j + 1])))
2186 if (j < npoints - 1)
2188 /* Found our position between points j and j+1.
2189 * Interpolate: x is the amount from j+1, (1-x) from point j
2190 * First treat possible jumps in temperature as a special case.
2192 if ((inputrec.opts.anneal_time[temperatureGroup][j + 1]
2193 - inputrec.opts.anneal_time[temperatureGroup][j])
2194 < GMX_REAL_EPS * 100)
2196 ir->opts.ref_t[temperatureGroup] = inputrec.opts.anneal_temp[temperatureGroup][j + 1];
2200 const real x = ((thist - inputrec.opts.anneal_time[temperatureGroup][j])
2201 / (inputrec.opts.anneal_time[temperatureGroup][j + 1]
2202 - inputrec.opts.anneal_time[temperatureGroup][j]));
2203 ir->opts.ref_t[temperatureGroup] =
2204 x * inputrec.opts.anneal_temp[temperatureGroup][j + 1]
2205 + (1 - x) * inputrec.opts.anneal_temp[temperatureGroup][j];
2210 ir->opts.ref_t[temperatureGroup] = inputrec.opts.anneal_temp[temperatureGroup][npoints - 1];
2214 upd->update_temperature_constants(*ir);
2217 void pleaseCiteCouplingAlgorithms(FILE* fplog, const t_inputrec& ir)
2219 if (EI_DYNAMICS(ir.eI))
2221 if (ir.etc == TemperatureCoupling::Berendsen)
2223 please_cite(fplog, "Berendsen84a");
2225 if (ir.etc == TemperatureCoupling::VRescale)
2227 please_cite(fplog, "Bussi2007a");
2229 if (ir.epc == PressureCoupling::CRescale)
2231 please_cite(fplog, "Bernetti2020");
2233 // TODO this is actually an integrator, not a coupling algorithm
2234 if (ir.eI == IntegrationAlgorithm::SD1)
2236 please_cite(fplog, "Goga2012");