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/mdatom.h"
65 #include "gromacs/mdtypes/state.h"
66 #include "gromacs/pbcutil/boxutilities.h"
67 #include "gromacs/pbcutil/pbc.h"
68 #include "gromacs/random/gammadistribution.h"
69 #include "gromacs/random/normaldistribution.h"
70 #include "gromacs/random/tabulatednormaldistribution.h"
71 #include "gromacs/random/threefry.h"
72 #include "gromacs/random/uniformrealdistribution.h"
73 #include "gromacs/utility/cstringutil.h"
74 #include "gromacs/utility/fatalerror.h"
75 #include "gromacs/utility/pleasecite.h"
76 #include "gromacs/utility/smalloc.h"
78 #define NTROTTERPARTS 3
80 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
82 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
83 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
85 #define MAX_SUZUKI_YOSHIDA_NUM 5
86 #define SUZUKI_YOSHIDA_NUM 5
88 static const double sy_const_1[] = { 1. };
89 static const double sy_const_3[] = { 0.828981543588751, -0.657963087177502, 0.828981543588751 };
90 static const double sy_const_5[] = { 0.2967324292201065,
96 static const double* sy_const[] = { nullptr, sy_const_1, nullptr, 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,
107 // This condition was explicitly checked in previous version, but should have never been satisfied
108 GMX_ASSERT(!(EI_VV(inputrec->eI)
109 && (inputrecNvtTrotter(inputrec) || inputrecNptTrotter(inputrec)
110 || inputrecNphTrotter(inputrec))),
111 "Temperature coupling was requested with velocity verlet and trotter");
113 bool doTemperatureCoupling = false;
115 // For VV temperature coupling parameters are updated on the current
116 // step, for the others - one step before.
117 if (inputrec->etc == etcNO)
119 doTemperatureCoupling = false;
121 else if (EI_VV(inputrec->eI))
123 doTemperatureCoupling = do_per_step(step, inputrec->nsttcouple);
127 doTemperatureCoupling = do_per_step(step + inputrec->nsttcouple - 1, inputrec->nsttcouple);
130 if (doTemperatureCoupling)
132 real dttc = inputrec->nsttcouple * inputrec->delta_t;
134 // TODO: berendsen_tcoupl(...), nosehoover_tcoupl(...) and vrescale_tcoupl(...) update
135 // temperature coupling parameters, which should be reflected in the name of these
137 switch (inputrec->etc)
141 berendsen_tcoupl(inputrec, ekind, dttc, state->therm_integral);
144 nosehoover_tcoupl(&(inputrec->opts),
147 state->nosehoover_xi.data(),
148 state->nosehoover_vxi.data(),
152 vrescale_tcoupl(inputrec, step, ekind, dttc, state->therm_integral.data());
155 /* rescale in place here */
156 if (EI_VV(inputrec->eI))
158 rescale_velocities(ekind, md, 0, md->homenr, state->v.rvec_array());
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 == epcPARRINELLORAHMAN
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 const matrix pressure,
198 const matrix forceVirial,
199 const matrix constraintVirial,
200 matrix pressureCouplingMu,
203 gmx::BoxDeformation* boxDeformation,
204 const bool scaleCoordinates)
207 int homenr = md->homenr;
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)
218 if (do_per_step(step, inputrec->nstpcouple))
220 real dtpc = inputrec->nstpcouple * dt;
221 pressureCouplingCalculateScalingMatrix<epcBERENDSEN>(fplog,
230 &state->baros_integral);
231 pressureCouplingScaleBoxAndCoordinates<epcBERENDSEN>(inputrec,
237 state->x.rvec_array(),
245 if (do_per_step(step, inputrec->nstpcouple))
247 real dtpc = inputrec->nstpcouple * dt;
248 pressureCouplingCalculateScalingMatrix<epcCRESCALE>(fplog,
257 &state->baros_integral);
258 pressureCouplingScaleBoxAndCoordinates<epcCRESCALE>(inputrec,
264 state->x.rvec_array(),
265 state->v.rvec_array(),
271 case (epcPARRINELLORAHMAN):
272 if (do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
274 /* The box velocities were updated in do_pr_pcoupl,
275 * but we dont change the box vectors until we get here
276 * since we need to be able to shift/unshift above.
278 real dtpc = inputrec->nstpcouple * dt;
279 for (int i = 0; i < DIM; i++)
281 for (int m = 0; m <= i; m++)
283 state->box[i][m] += dtpc * state->boxv[i][m];
286 preserve_box_shape(inputrec, state->box_rel, state->box);
288 /* Scale the coordinates */
289 if (scaleCoordinates)
291 auto x = state->x.rvec_array();
292 for (int n = start; n < start + homenr; n++)
294 tmvmul_ur0(pressureCouplingMu, x[n], x[n]);
300 switch (inputrec->epct)
302 case (epctISOTROPIC):
303 /* DIM * eta = ln V. so DIM*eta_new = DIM*eta_old + DIM*dt*veta =>
304 ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) =>
305 Side length scales as exp(veta*dt) */
307 msmul(state->box, std::exp(state->veta * dt), state->box);
309 /* Relate veta to boxv. veta = d(eta)/dT = (1/DIM)*1/V dV/dT.
310 o If we assume isotropic scaling, and box length scaling
311 factor L, then V = L^DIM (det(M)). So dV/dt = DIM
312 L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt. The
313 determinant of B is L^DIM det(M), and the determinant
314 of dB/dt is (dL/dT)^DIM det (M). veta will be
315 (det(dB/dT)/det(B))^(1/3). Then since M =
316 B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */
318 msmul(state->box, state->veta, state->boxv);
328 auto localX = makeArrayRef(state->x).subArray(start, homenr);
329 boxDeformation->apply(localX, state->box, step);
333 extern gmx_bool update_randomize_velocities(const t_inputrec* ir,
337 gmx::ArrayRef<gmx::RVec> v,
338 const gmx::Update* upd,
339 const gmx::Constraints* constr)
342 real rate = (ir->delta_t) / ir->opts.tau_t[0];
344 if (ir->etc == etcANDERSEN && constr != nullptr)
346 /* Currently, Andersen thermostat does not support constrained
347 systems. Functionality exists in the andersen_tcoupl
348 function in GROMACS 4.5.7 to allow this combination. That
349 code could be ported to the current random-number
350 generation approach, but has not yet been done because of
351 lack of time and resources. */
353 "Normal Andersen is currently not supported with constraints, use massive "
357 /* proceed with andersen if 1) it's fixed probability per
358 particle andersen or 2) it's massive andersen and it's tau_t/dt */
359 if ((ir->etc == etcANDERSEN) || do_per_step(step, gmx::roundToInt(1.0 / rate)))
362 ir, step, cr, md, v, rate, upd->getAndersenRandomizeGroup(), upd->getBoltzmanFactor());
369 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
373 {0.828981543588751,-0.657963087177502,0.828981543588751},
375 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
378 /* these integration routines are only referenced inside this file */
379 static void NHC_trotter(const t_grpopts* opts,
381 const gmx_ekindata_t* ekind,
387 const t_extmass* MassQ,
388 gmx_bool bEkinAveVel)
391 /* general routine for both barostat and thermostat nose hoover chains */
394 double Ekin, Efac, reft, kT, nd;
399 int mstepsi, mstepsj;
400 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
401 int nh = opts->nhchainlength;
404 mstepsi = mstepsj = ns;
406 /* if scalefac is NULL, we are doing the NHC of the barostat */
409 if (scalefac == nullptr)
414 for (i = 0; i < nvar; i++)
417 /* make it easier to iterate by selecting
418 out the sub-array that corresponds to this T group */
422 gmx::ArrayRef<const double> iQinv;
425 iQinv = gmx::arrayRefFromArray(&MassQ->QPinv[i * nh], nh);
426 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
427 reft = std::max<real>(0, opts->ref_t[0]);
428 Ekin = gmx::square(*veta) / MassQ->Winv;
432 iQinv = gmx::arrayRefFromArray(&MassQ->Qinv[i * nh], nh);
433 const t_grp_tcstat* tcstat = &ekind->tcstat[i];
435 reft = std::max<real>(0, opts->ref_t[i]);
438 Ekin = 2 * trace(tcstat->ekinf) * tcstat->ekinscalef_nhc;
442 Ekin = 2 * trace(tcstat->ekinh) * tcstat->ekinscaleh_nhc;
447 for (mi = 0; mi < mstepsi; mi++)
449 for (mj = 0; mj < mstepsj; mj++)
451 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
452 dt = sy_const[ns][mj] * dtfull / mstepsi;
454 /* compute the thermal forces */
455 GQ[0] = iQinv[0] * (Ekin - nd * kT);
457 for (j = 0; j < nh - 1; j++)
459 if (iQinv[j + 1] > 0)
461 /* we actually don't need to update here if we save the
462 state of the GQ, but it's easier to just recompute*/
463 GQ[j + 1] = iQinv[j + 1] * ((gmx::square(ivxi[j]) / iQinv[j]) - kT);
471 ivxi[nh - 1] += 0.25 * dt * GQ[nh - 1];
472 for (j = nh - 1; j > 0; j--)
474 Efac = exp(-0.125 * dt * ivxi[j]);
475 ivxi[j - 1] = Efac * (ivxi[j - 1] * Efac + 0.25 * dt * GQ[j - 1]);
478 Efac = exp(-0.5 * dt * ivxi[0]);
487 Ekin *= (Efac * Efac);
489 /* Issue - if the KE is an average of the last and the current temperatures, then we
490 might not be able to scale the kinetic energy directly with this factor. Might
491 take more bookkeeping -- have to think about this a bit more . . . */
493 GQ[0] = iQinv[0] * (Ekin - nd * kT);
495 /* update thermostat positions */
496 for (j = 0; j < nh; j++)
498 ixi[j] += 0.5 * dt * ivxi[j];
501 for (j = 0; j < nh - 1; j++)
503 Efac = exp(-0.125 * dt * ivxi[j + 1]);
504 ivxi[j] = Efac * (ivxi[j] * Efac + 0.25 * dt * GQ[j]);
505 if (iQinv[j + 1] > 0)
507 GQ[j + 1] = iQinv[j + 1] * ((gmx::square(ivxi[j]) / iQinv[j]) - kT);
514 ivxi[nh - 1] += 0.25 * dt * GQ[nh - 1];
521 static void boxv_trotter(const t_inputrec* ir,
525 const gmx_ekindata_t* ekind,
528 const t_extmass* MassQ)
535 tensor ekinmod, localpres;
537 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
538 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
541 if (ir->epct == epctSEMIISOTROPIC)
550 /* eta is in pure units. veta is in units of ps^-1. GW is in
551 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
552 taken to use only RATIOS of eta in updating the volume. */
554 /* we take the partial pressure tensors, modify the
555 kinetic energy tensor, and recovert to pressure */
557 if (ir->opts.nrdf[0] == 0)
559 gmx_fatal(FARGS, "Barostat is coupled to a T-group with no degrees of freedom\n");
561 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
562 alpha = 1.0 + DIM / (static_cast<double>(ir->opts.nrdf[0]));
563 alpha *= ekind->tcstat[0].ekinscalef_nhc;
564 msmul(ekind->ekin, alpha, ekinmod);
565 /* for now, we use Elr = 0, because if you want to get it right, you
566 really should be using PME. Maybe print a warning? */
568 pscal = calc_pres(ir->pbcType, nwall, box, ekinmod, vir, localpres) + pcorr;
571 GW = (vol * (MassQ->Winv / PRESFAC)) * (DIM * pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
573 *veta += 0.5 * dt * GW;
577 * This file implements temperature and pressure coupling algorithms:
578 * For now only the Weak coupling and the modified weak coupling.
580 * Furthermore computation of pressure and temperature is done here
584 real calc_pres(PbcType pbcType, int nwall, const matrix box, const tensor ekin, const tensor vir, tensor pres)
589 if (pbcType == PbcType::No || (pbcType == PbcType::XY && nwall != 2))
595 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
596 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
600 fac = PRESFAC * 2.0 / det(box);
601 for (n = 0; (n < DIM); n++)
603 for (m = 0; (m < DIM); m++)
605 pres[n][m] = (ekin[n][m] - vir[n][m]) * fac;
611 pr_rvecs(debug, 0, "PC: pres", pres, DIM);
612 pr_rvecs(debug, 0, "PC: ekin", ekin, DIM);
613 pr_rvecs(debug, 0, "PC: vir ", vir, DIM);
614 pr_rvecs(debug, 0, "PC: box ", box, DIM);
617 return trace(pres) / DIM;
620 real calc_temp(real ekin, real nrdf)
624 return (2.0 * ekin) / (nrdf * BOLTZ);
632 /*! \brief Sets 1/mass for Parrinello-Rahman in wInv; NOTE: PRESFAC is not included, so not in GROMACS units! */
633 static void calcParrinelloRahmanInvMass(const t_inputrec* ir, const matrix box, tensor wInv)
637 /* TODO: See if we can make the mass independent of the box size */
638 maxBoxLength = std::max(box[XX][XX], box[YY][YY]);
639 maxBoxLength = std::max(maxBoxLength, box[ZZ][ZZ]);
641 for (int d = 0; d < DIM; d++)
643 for (int n = 0; n < DIM; n++)
645 wInv[d][n] = (4 * M_PI * M_PI * ir->compress[d][n]) / (3 * ir->tau_p * ir->tau_p * maxBoxLength);
650 void parrinellorahman_pcoupl(FILE* fplog,
652 const t_inputrec* ir,
662 /* This doesn't do any coordinate updating. It just
663 * integrates the box vector equations from the calculated
664 * acceleration due to pressure difference. We also compute
665 * the tensor M which is used in update to couple the particle
666 * coordinates to the box vectors.
668 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
671 * M_nk = (h') * (h' * h + h' h) * h
673 * with the dots denoting time derivatives and h is the transformation from
674 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
675 * This also goes for the pressure and M tensors - they are transposed relative
676 * to ours. Our equation thus becomes:
679 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
681 * where b is the gromacs box matrix.
682 * Our box accelerations are given by
684 * b = vol/W inv(box') * (P-ref_P) (=h')
687 real vol = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
688 real atot, arel, change;
689 tensor invbox, pdiff, t1, t2;
691 gmx::invertBoxMatrix(box, invbox);
695 /* Note that PRESFAC does not occur here.
696 * The pressure and compressibility always occur as a product,
697 * therefore the pressure unit drops out.
700 calcParrinelloRahmanInvMass(ir, box, winv);
702 m_sub(pres, ir->ref_p, pdiff);
704 if (ir->epct == epctSURFACETENSION)
706 /* Unlike Berendsen coupling it might not be trivial to include a z
707 * pressure correction here? On the other hand we don't scale the
708 * box momentarily, but change accelerations, so it might not be crucial.
710 real xy_pressure = 0.5 * (pres[XX][XX] + pres[YY][YY]);
711 for (int d = 0; d < ZZ; d++)
713 pdiff[d][d] = (xy_pressure - (pres[ZZ][ZZ] - ir->ref_p[d][d] / box[d][d]));
717 tmmul(invbox, pdiff, t1);
718 /* Move the off-diagonal elements of the 'force' to one side to ensure
719 * that we obey the box constraints.
721 for (int d = 0; d < DIM; d++)
723 for (int n = 0; n < d; n++)
725 t1[d][n] += t1[n][d];
732 case epctANISOTROPIC:
733 for (int d = 0; d < DIM; d++)
735 for (int n = 0; n <= d; n++)
737 t1[d][n] *= winv[d][n] * vol;
742 /* calculate total volume acceleration */
743 atot = box[XX][XX] * box[YY][YY] * t1[ZZ][ZZ] + box[XX][XX] * t1[YY][YY] * box[ZZ][ZZ]
744 + t1[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
745 arel = atot / (3 * vol);
746 /* set all RELATIVE box accelerations equal, and maintain total V
748 for (int d = 0; d < DIM; d++)
750 for (int n = 0; n <= d; n++)
752 t1[d][n] = winv[0][0] * vol * arel * box[d][n];
756 case epctSEMIISOTROPIC:
757 case epctSURFACETENSION:
758 /* Note the correction to pdiff above for surftens. coupling */
760 /* calculate total XY volume acceleration */
761 atot = box[XX][XX] * t1[YY][YY] + t1[XX][XX] * box[YY][YY];
762 arel = atot / (2 * box[XX][XX] * box[YY][YY]);
763 /* set RELATIVE XY box accelerations equal, and maintain total V
764 * change speed. Dont change the third box vector accelerations */
765 for (int d = 0; d < ZZ; d++)
767 for (int n = 0; n <= d; n++)
769 t1[d][n] = winv[d][n] * vol * arel * box[d][n];
772 for (int n = 0; n < DIM; n++)
774 t1[ZZ][n] *= winv[ZZ][n] * vol;
779 "Parrinello-Rahman pressure coupling type %s "
780 "not supported yet\n",
781 EPCOUPLTYPETYPE(ir->epct));
785 for (int d = 0; d < DIM; d++)
787 for (int n = 0; n <= d; n++)
789 boxv[d][n] += dt * t1[d][n];
791 /* We do NOT update the box vectors themselves here, since
792 * we need them for shifting later. It is instead done last
793 * in the update() routine.
796 /* Calculate the change relative to diagonal elements-
797 since it's perfectly ok for the off-diagonal ones to
798 be zero it doesn't make sense to check the change relative
802 change = std::fabs(dt * boxv[d][n] / box[d][d]);
804 if (change > maxchange)
811 if (maxchange > 0.01 && fplog)
815 "\nStep %s Warning: Pressure scaling more than 1%%. "
816 "This may mean your system\n is not yet equilibrated. "
817 "Use of Parrinello-Rahman pressure coupling during\n"
818 "equilibration can lead to simulation instability, "
819 "and is discouraged.\n",
820 gmx_step_str(step, buf));
824 preserve_box_shape(ir, box_rel, boxv);
826 mtmul(boxv, box, t1); /* t1=boxv * b' */
827 mmul(invbox, t1, t2);
828 mtmul(t2, invbox, M);
830 /* Determine the scaling matrix mu for the coordinates */
831 for (int d = 0; d < DIM; d++)
833 for (int n = 0; n <= d; n++)
835 t1[d][n] = box[d][n] + dt * boxv[d][n];
838 preserve_box_shape(ir, box_rel, t1);
839 /* t1 is the box at t+dt, determine mu as the relative change */
840 mmul_ur0(invbox, t1, mu);
843 //! Return compressibility factor for entry (i,j) of Berendsen / C-rescale scaling matrix
844 static inline real compressibilityFactor(int i, int j, const t_inputrec* ir, real dt)
846 return ir->compress[i][j] * dt / ir->tau_p;
849 //! Details of Berendsen / C-rescale scaling matrix calculation
850 template<int pressureCouplingType>
851 static void calculateScalingMatrixImplDetail(const t_inputrec* ir,
856 real scalar_pressure,
860 //! Calculate Berendsen / C-rescale scaling matrix
861 template<int pressureCouplingType>
862 static void calculateScalingMatrixImpl(const t_inputrec* ir,
869 real scalar_pressure = 0;
870 real xy_pressure = 0;
871 for (int d = 0; d < DIM; d++)
873 scalar_pressure += pres[d][d] / DIM;
876 xy_pressure += pres[d][d] / (DIM - 1);
880 calculateScalingMatrixImplDetail<pressureCouplingType>(
881 ir, mu, dt, pres, box, scalar_pressure, xy_pressure, step);
885 void calculateScalingMatrixImplDetail<epcBERENDSEN>(const t_inputrec* ir,
890 real scalar_pressure,
892 int64_t gmx_unused step)
898 for (int d = 0; d < DIM; d++)
900 mu[d][d] = 1.0 - compressibilityFactor(d, d, ir, dt) * (ir->ref_p[d][d] - scalar_pressure) / DIM;
903 case epctSEMIISOTROPIC:
904 for (int d = 0; d < ZZ; d++)
906 mu[d][d] = 1.0 - compressibilityFactor(d, d, ir, dt) * (ir->ref_p[d][d] - xy_pressure) / DIM;
909 1.0 - compressibilityFactor(ZZ, ZZ, ir, dt) * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM;
911 case epctANISOTROPIC:
912 for (int d = 0; d < DIM; d++)
914 for (int n = 0; n < DIM; n++)
916 mu[d][n] = (d == n ? 1.0 : 0.0)
917 - compressibilityFactor(d, n, ir, dt) * (ir->ref_p[d][n] - pres[d][n]) / DIM;
921 case epctSURFACETENSION:
922 /* ir->ref_p[0/1] is the reference surface-tension times *
923 * the number of surfaces */
924 if (ir->compress[ZZ][ZZ] != 0.0F)
926 p_corr_z = dt / ir->tau_p * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
930 /* when the compressibity is zero, set the pressure correction *
931 * in the z-direction to zero to get the correct surface tension */
934 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ] * p_corr_z;
935 for (int d = 0; d < DIM - 1; d++)
938 + compressibilityFactor(d, d, ir, dt)
939 * (ir->ref_p[d][d] / (mu[ZZ][ZZ] * box[ZZ][ZZ])
940 - (pres[ZZ][ZZ] + p_corr_z - xy_pressure))
946 "Berendsen pressure coupling type %s not supported yet\n",
947 EPCOUPLTYPETYPE(ir->epct));
952 void calculateScalingMatrixImplDetail<epcCRESCALE>(const t_inputrec* ir,
957 real scalar_pressure,
961 gmx::ThreeFry2x64<64> rng(ir->ld_seed, gmx::RandomDomain::Barostat);
962 gmx::NormalDistribution<real> normalDist;
963 rng.restart(step, 0);
965 for (int d = 0; d < DIM; d++)
971 real kt = ir->opts.ref_t[0] * BOLTZ;
980 gauss = normalDist(rng);
981 for (int d = 0; d < DIM; d++)
983 const real factor = compressibilityFactor(d, d, ir, dt);
984 mu[d][d] = std::exp(-factor * (ir->ref_p[d][d] - scalar_pressure) / DIM
985 + std::sqrt(2.0 * kt * factor * PRESFAC / vol) * gauss / DIM);
988 case epctSEMIISOTROPIC:
989 gauss = normalDist(rng);
990 gauss2 = normalDist(rng);
991 for (int d = 0; d < ZZ; d++)
993 const real factor = compressibilityFactor(d, d, ir, dt);
994 mu[d][d] = std::exp(-factor * (ir->ref_p[d][d] - xy_pressure) / DIM
995 + std::sqrt((DIM - 1) * 2.0 * kt * factor * PRESFAC / vol / DIM)
996 / (DIM - 1) * gauss);
999 const real factor = compressibilityFactor(ZZ, ZZ, ir, dt);
1000 mu[ZZ][ZZ] = std::exp(-factor * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM
1001 + std::sqrt(2.0 * kt * factor * PRESFAC / vol / DIM) * gauss2);
1004 case epctSURFACETENSION:
1005 gauss = normalDist(rng);
1006 gauss2 = normalDist(rng);
1007 for (int d = 0; d < ZZ; d++)
1009 const real factor = compressibilityFactor(d, d, ir, dt);
1010 /* Notice: we here use ref_p[ZZ][ZZ] as isotropic pressure and ir->ref_p[d][d] as surface tension */
1011 mu[d][d] = std::exp(
1012 -factor * (ir->ref_p[ZZ][ZZ] - ir->ref_p[d][d] / box[ZZ][ZZ] - xy_pressure) / DIM
1013 + std::sqrt(4.0 / 3.0 * kt * factor * PRESFAC / vol) / (DIM - 1) * gauss);
1016 const real factor = compressibilityFactor(ZZ, ZZ, ir, dt);
1017 mu[ZZ][ZZ] = std::exp(-factor * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM
1018 + std::sqrt(2.0 / 3.0 * kt * factor * PRESFAC / vol) * gauss2);
1023 "C-rescale pressure coupling type %s not supported yet\n",
1024 EPCOUPLTYPETYPE(ir->epct));
1028 template<int pressureCouplingType>
1029 void pressureCouplingCalculateScalingMatrix(FILE* fplog,
1031 const t_inputrec* ir,
1035 const matrix force_vir,
1036 const matrix constraint_vir,
1038 double* baros_integral)
1040 static_assert(pressureCouplingType == epcBERENDSEN || pressureCouplingType == epcCRESCALE,
1041 "pressureCouplingCalculateScalingMatrix is only implemented for Berendsen and "
1042 "C-rescale pressure coupling");
1044 calculateScalingMatrixImpl<pressureCouplingType>(ir, mu, dt, pres, box, step);
1046 /* To fulfill the orientation restrictions on triclinic boxes
1047 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
1048 * the other elements of mu to first order.
1050 mu[YY][XX] += mu[XX][YY];
1051 mu[ZZ][XX] += mu[XX][ZZ];
1052 mu[ZZ][YY] += mu[YY][ZZ];
1057 /* Keep track of the work the barostat applies on the system.
1058 * Without constraints force_vir tells us how Epot changes when scaling.
1059 * With constraints constraint_vir gives us the constraint contribution
1060 * to both Epot and Ekin. Although we are not scaling velocities, scaling
1061 * the coordinates leads to scaling of distances involved in constraints.
1062 * This in turn changes the angular momentum (even if the constrained
1063 * distances are corrected at the next step). The kinetic component
1064 * of the constraint virial captures the angular momentum change.
1066 for (int d = 0; d < DIM; d++)
1068 for (int n = 0; n <= d; n++)
1071 2 * (mu[d][n] - (n == d ? 1 : 0)) * (force_vir[d][n] + constraint_vir[d][n]);
1077 pr_rvecs(debug, 0, "PC: pres ", pres, 3);
1078 pr_rvecs(debug, 0, "PC: mu ", mu, 3);
1081 if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 || mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01
1082 || mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
1087 "\nStep %s Warning: pressure scaling more than 1%%, "
1089 gmx_step_str(step, buf2),
1095 fprintf(fplog, "%s", buf);
1097 fprintf(stderr, "%s", buf);
1101 template<int pressureCouplingType>
1102 void pressureCouplingScaleBoxAndCoordinates(const t_inputrec* ir,
1110 const unsigned short cFREEZE[],
1112 const bool scaleCoordinates)
1114 static_assert(pressureCouplingType == epcBERENDSEN || pressureCouplingType == epcCRESCALE,
1115 "pressureCouplingScaleBoxAndCoordinates is only implemented for Berendsen and "
1116 "C-rescale pressure coupling");
1118 ivec* nFreeze = ir->opts.nFreeze;
1120 if (pressureCouplingType == epcCRESCALE)
1122 gmx::invertBoxMatrix(mu, inv_mu);
1125 /* Scale the positions and the velocities */
1126 if (scaleCoordinates)
1128 const int gmx_unused numThreads = gmx_omp_nthreads_get(emntUpdate);
1129 #pragma omp parallel for num_threads(numThreads) schedule(static)
1130 for (int n = start; n < start + nr_atoms; n++)
1132 // Trivial OpenMP region that does not throw
1134 if (cFREEZE != nullptr)
1139 if (!nFreeze[g][XX])
1141 x[n][XX] = mu[XX][XX] * x[n][XX] + mu[YY][XX] * x[n][YY] + mu[ZZ][XX] * x[n][ZZ];
1142 if (pressureCouplingType == epcCRESCALE)
1144 v[n][XX] = inv_mu[XX][XX] * v[n][XX] + inv_mu[YY][XX] * v[n][YY]
1145 + inv_mu[ZZ][XX] * v[n][ZZ];
1148 if (!nFreeze[g][YY])
1150 x[n][YY] = mu[YY][YY] * x[n][YY] + mu[ZZ][YY] * x[n][ZZ];
1151 if (pressureCouplingType == epcCRESCALE)
1153 v[n][YY] = inv_mu[YY][YY] * v[n][YY] + inv_mu[ZZ][YY] * v[n][ZZ];
1156 if (!nFreeze[g][ZZ])
1158 x[n][ZZ] = mu[ZZ][ZZ] * x[n][ZZ];
1159 if (pressureCouplingType == epcCRESCALE)
1161 v[n][ZZ] = inv_mu[ZZ][ZZ] * v[n][ZZ];
1166 /* compute final boxlengths */
1167 for (int d = 0; d < DIM; d++)
1169 box[d][XX] = mu[XX][XX] * box[d][XX] + mu[YY][XX] * box[d][YY] + mu[ZZ][XX] * box[d][ZZ];
1170 box[d][YY] = mu[YY][YY] * box[d][YY] + mu[ZZ][YY] * box[d][ZZ];
1171 box[d][ZZ] = mu[ZZ][ZZ] * box[d][ZZ];
1174 preserve_box_shape(ir, box_rel, box);
1176 /* (un)shifting should NOT be done after this,
1177 * since the box vectors might have changed
1179 inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
1182 void berendsen_tcoupl(const t_inputrec* ir, gmx_ekindata_t* ekind, real dt, std::vector<double>& therm_integral)
1184 const t_grpopts* opts = &ir->opts;
1186 for (int i = 0; (i < opts->ngtc); i++)
1192 Ek = trace(ekind->tcstat[i].ekinf);
1193 T = ekind->tcstat[i].T;
1197 Ek = trace(ekind->tcstat[i].ekinh);
1198 T = ekind->tcstat[i].Th;
1201 if ((opts->tau_t[i] > 0) && (T > 0.0))
1203 real reft = std::max<real>(0, opts->ref_t[i]);
1204 real lll = std::sqrt(1.0 + (dt / opts->tau_t[i]) * (reft / T - 1.0));
1205 ekind->tcstat[i].lambda = std::max<real>(std::min<real>(lll, 1.25), 0.8);
1209 ekind->tcstat[i].lambda = 1.0;
1212 /* Keep track of the amount of energy we are adding to the system */
1213 therm_integral[i] -= (gmx::square(ekind->tcstat[i].lambda) - 1) * Ek;
1217 fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n", i, T, ekind->tcstat[i].lambda);
1222 void andersen_tcoupl(const t_inputrec* ir,
1224 const t_commrec* cr,
1225 const t_mdatoms* md,
1226 gmx::ArrayRef<gmx::RVec> v,
1228 const std::vector<bool>& randomize,
1229 gmx::ArrayRef<const real> boltzfac)
1231 const int* gatindex = (DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1234 gmx::ThreeFry2x64<0> rng(ir->andersen_seed, gmx::RandomDomain::Thermostat);
1235 gmx::UniformRealDistribution<real> uniformDist;
1236 gmx::TabulatedNormalDistribution<real, 14> normalDist;
1238 /* randomize the velocities of the selected particles */
1240 for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */
1242 int ng = gatindex ? gatindex[i] : i;
1243 gmx_bool bRandomize;
1245 rng.restart(step, ng);
1249 gc = md->cTC[i]; /* assign the atom to a temperature group if there are more than one */
1253 if (ir->etc == etcANDERSENMASSIVE)
1255 /* Randomize particle always */
1260 /* Randomize particle probabilistically */
1261 uniformDist.reset();
1262 bRandomize = uniformDist(rng) < rate;
1269 scal = std::sqrt(boltzfac[gc] * md->invmass[i]);
1273 for (d = 0; d < DIM; d++)
1275 v[i][d] = scal * normalDist(rng);
1283 void nosehoover_tcoupl(const t_grpopts* opts,
1284 const gmx_ekindata_t* ekind,
1288 const t_extmass* MassQ)
1293 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
1295 for (i = 0; (i < opts->ngtc); i++)
1297 reft = std::max<real>(0, opts->ref_t[i]);
1299 vxi[i] += dt * MassQ->Qinv[i] * (ekind->tcstat[i].Th - reft);
1300 xi[i] += dt * (oldvxi + vxi[i]) * 0.5;
1304 void trotter_update(const t_inputrec* ir,
1306 gmx_ekindata_t* ekind,
1307 const gmx_enerdata_t* enerd,
1310 const t_mdatoms* md,
1311 const t_extmass* MassQ,
1312 gmx::ArrayRef<std::vector<int>> trotter_seqlist,
1316 int n, i, d, ngtc, gc = 0, t;
1317 t_grp_tcstat* tcstat;
1318 const t_grpopts* opts;
1321 double * scalefac, dtc;
1322 rvec sumv = { 0, 0, 0 };
1325 if (trotter_seqno <= ettTSEQ2)
1327 step_eff = step - 1; /* the velocity verlet calls are actually out of order -- the first
1328 half step is actually the last half step from the previous step.
1329 Thus the first half step actually corresponds to the n-1 step*/
1336 bCouple = (ir->nsttcouple == 1 || do_per_step(step_eff + ir->nsttcouple, ir->nsttcouple));
1338 const gmx::ArrayRef<const int> trotter_seq = trotter_seqlist[trotter_seqno];
1340 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
1344 dtc = ir->nsttcouple * ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
1345 opts = &(ir->opts); /* just for ease of referencing */
1348 snew(scalefac, opts->ngtc);
1349 for (i = 0; i < ngtc; i++)
1353 /* execute the series of trotter updates specified in the trotterpart array */
1355 for (i = 0; i < NTROTTERPARTS; i++)
1357 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
1358 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2)
1359 || (trotter_seq[i] == etrtNHC2))
1368 auto v = makeArrayRef(state->v);
1369 switch (trotter_seq[i])
1373 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir, enerd->term[F_PDISPCORR], MassQ);
1381 state->nhpres_xi.data(),
1382 state->nhpres_vxi.data(),
1394 state->nosehoover_xi.data(),
1395 state->nosehoover_vxi.data(),
1400 /* need to rescale the kinetic energies and velocities here. Could
1401 scale the velocities later, but we need them scaled in order to
1402 produce the correct outputs, so we'll scale them here. */
1404 for (t = 0; t < ngtc; t++)
1406 tcstat = &ekind->tcstat[t];
1407 tcstat->vscale_nhc = scalefac[t];
1408 tcstat->ekinscaleh_nhc *= (scalefac[t] * scalefac[t]);
1409 tcstat->ekinscalef_nhc *= (scalefac[t] * scalefac[t]);
1411 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
1412 /* but do we actually need the total? */
1414 /* modify the velocities as well */
1415 for (n = 0; n < md->homenr; n++)
1417 if (md->cTC) /* does this conditional need to be here? is this always true?*/
1421 for (d = 0; d < DIM; d++)
1423 v[n][d] *= scalefac[gc];
1428 for (d = 0; d < DIM; d++)
1430 sumv[d] += (v[n][d]) / md->invmass[n];
1438 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
1443 extern void init_npt_masses(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool bInit)
1445 int n, i, j, d, ngtc, nh;
1446 const t_grpopts* opts;
1447 real reft, kT, ndj, nd;
1449 opts = &(ir->opts); /* just for ease of referencing */
1450 ngtc = ir->opts.ngtc;
1451 nh = state->nhchainlength;
1457 MassQ->Qinv.resize(ngtc);
1459 for (i = 0; (i < ngtc); i++)
1461 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1463 MassQ->Qinv[i] = 1.0 / (gmx::square(opts->tau_t[i] / M_2PI) * opts->ref_t[i]);
1467 MassQ->Qinv[i] = 0.0;
1471 else if (EI_VV(ir->eI))
1473 /* Set pressure variables */
1477 if (state->vol0 == 0)
1479 state->vol0 = det(state->box);
1480 /* because we start by defining a fixed
1481 compressibility, we need the volume at this
1482 compressibility to solve the problem. */
1486 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1487 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1488 MassQ->Winv = (PRESFAC * trace(ir->compress) * BOLTZ * opts->ref_t[0])
1489 / (DIM * state->vol0 * gmx::square(ir->tau_p / M_2PI));
1490 /* An alternate mass definition, from Tuckerman et al. */
1491 /* MassQ->Winv = 1.0/(gmx::square(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1492 for (d = 0; d < DIM; d++)
1494 for (n = 0; n < DIM; n++)
1496 MassQ->Winvm[d][n] =
1497 PRESFAC * ir->compress[d][n] / (state->vol0 * gmx::square(ir->tau_p / M_2PI));
1498 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1499 before using MTTK for anisotropic states.*/
1502 /* Allocate space for thermostat variables */
1505 MassQ->Qinv.resize(ngtc * nh);
1508 /* now, set temperature variables */
1509 for (i = 0; i < ngtc; i++)
1511 if (opts->tau_t[i] > 0 && opts->ref_t[i] > 0 && opts->nrdf[i] > 0)
1513 reft = std::max<real>(0, opts->ref_t[i]);
1516 for (j = 0; j < nh; j++)
1526 MassQ->Qinv[i * nh + j] = 1.0 / (gmx::square(opts->tau_t[i] / M_2PI) * ndj * kT);
1531 for (j = 0; j < nh; j++)
1533 MassQ->Qinv[i * nh + j] = 0.0;
1540 std::array<std::vector<int>, ettTSEQMAX>
1541 init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool bTrotter)
1543 int i, j, nnhpres, nh;
1544 const t_grpopts* opts;
1545 real bmass, qmass, reft, kT;
1547 opts = &(ir->opts); /* just for ease of referencing */
1548 nnhpres = state->nnhpres;
1549 nh = state->nhchainlength;
1551 if (EI_VV(ir->eI) && (ir->epc == epcMTTK) && (ir->etc != etcNOSEHOOVER))
1553 gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
1556 init_npt_masses(ir, state, MassQ, TRUE);
1558 /* first, initialize clear all the trotter calls */
1559 std::array<std::vector<int>, ettTSEQMAX> trotter_seq;
1560 for (i = 0; i < ettTSEQMAX; i++)
1562 trotter_seq[i].resize(NTROTTERPARTS, etrtNONE);
1563 trotter_seq[i][0] = etrtSKIPALL;
1568 /* no trotter calls, so we never use the values in the array.
1569 * We access them (so we need to define them, but ignore
1575 /* compute the kinetic energy by using the half step velocities or
1576 * the kinetic energies, depending on the order of the trotter calls */
1580 if (inputrecNptTrotter(ir))
1582 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1583 We start with the initial one. */
1584 /* first, a round that estimates veta. */
1585 trotter_seq[0][0] = etrtBAROV;
1587 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1589 /* The first half trotter update */
1590 trotter_seq[2][0] = etrtBAROV;
1591 trotter_seq[2][1] = etrtNHC;
1592 trotter_seq[2][2] = etrtBARONHC;
1594 /* The second half trotter update */
1595 trotter_seq[3][0] = etrtBARONHC;
1596 trotter_seq[3][1] = etrtNHC;
1597 trotter_seq[3][2] = etrtBAROV;
1599 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1601 else if (inputrecNvtTrotter(ir))
1603 /* This is the easy version - there are only two calls, both the same.
1604 Otherwise, even easier -- no calls */
1605 trotter_seq[2][0] = etrtNHC;
1606 trotter_seq[3][0] = etrtNHC;
1608 else if (inputrecNphTrotter(ir))
1610 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1611 We start with the initial one. */
1612 /* first, a round that estimates veta. */
1613 trotter_seq[0][0] = etrtBAROV;
1615 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1617 /* The first half trotter update */
1618 trotter_seq[2][0] = etrtBAROV;
1619 trotter_seq[2][1] = etrtBARONHC;
1621 /* The second half trotter update */
1622 trotter_seq[3][0] = etrtBARONHC;
1623 trotter_seq[3][1] = etrtBAROV;
1625 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1628 else if (ir->eI == eiVVAK)
1630 if (inputrecNptTrotter(ir))
1632 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1633 We start with the initial one. */
1634 /* first, a round that estimates veta. */
1635 trotter_seq[0][0] = etrtBAROV;
1637 /* The first half trotter update, part 1 -- double update, because it commutes */
1638 trotter_seq[1][0] = etrtNHC;
1640 /* The first half trotter update, part 2 */
1641 trotter_seq[2][0] = etrtBAROV;
1642 trotter_seq[2][1] = etrtBARONHC;
1644 /* The second half trotter update, part 1 */
1645 trotter_seq[3][0] = etrtBARONHC;
1646 trotter_seq[3][1] = etrtBAROV;
1648 /* The second half trotter update */
1649 trotter_seq[4][0] = etrtNHC;
1651 else if (inputrecNvtTrotter(ir))
1653 /* This is the easy version - there is only one call, both the same.
1654 Otherwise, even easier -- no calls */
1655 trotter_seq[1][0] = etrtNHC;
1656 trotter_seq[4][0] = etrtNHC;
1658 else if (inputrecNphTrotter(ir))
1660 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1661 We start with the initial one. */
1662 /* first, a round that estimates veta. */
1663 trotter_seq[0][0] = etrtBAROV;
1665 /* The first half trotter update, part 1 -- leave zero */
1666 trotter_seq[1][0] = etrtNHC;
1668 /* The first half trotter update, part 2 */
1669 trotter_seq[2][0] = etrtBAROV;
1670 trotter_seq[2][1] = etrtBARONHC;
1672 /* The second half trotter update, part 1 */
1673 trotter_seq[3][0] = etrtBARONHC;
1674 trotter_seq[3][1] = etrtBAROV;
1676 /* The second half trotter update -- blank for now */
1683 default: bmass = DIM * DIM; /* recommended mass parameters for isotropic barostat */
1686 MassQ->QPinv.resize(nnhpres * opts->nhchainlength);
1688 /* barostat temperature */
1689 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1691 reft = std::max<real>(0, opts->ref_t[0]);
1693 for (i = 0; i < nnhpres; i++)
1695 for (j = 0; j < nh; j++)
1705 MassQ->QPinv[i * opts->nhchainlength + j] =
1706 1.0 / (gmx::square(opts->tau_t[0] / M_2PI) * qmass * kT);
1712 for (i = 0; i < nnhpres; i++)
1714 for (j = 0; j < nh; j++)
1716 MassQ->QPinv[i * nh + j] = 0.0;
1723 static real energyNoseHoover(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1727 int nh = state->nhchainlength;
1729 for (int i = 0; i < ir->opts.ngtc; i++)
1731 const double* ixi = &state->nosehoover_xi[i * nh];
1732 const double* ivxi = &state->nosehoover_vxi[i * nh];
1733 const double* iQinv = &(MassQ->Qinv[i * nh]);
1735 real nd = ir->opts.nrdf[i];
1736 real reft = std::max<real>(ir->opts.ref_t[i], 0);
1737 real kT = BOLTZ * reft;
1741 if (inputrecNvtTrotter(ir) || inputrecNptTrotter(ir))
1743 /* contribution from the thermal momenta of the NH chain */
1744 for (int j = 0; j < nh; j++)
1748 energy += 0.5 * gmx::square(ivxi[j]) / iQinv[j];
1749 /* contribution from the thermal variable of the NH chain */
1759 energy += ndj * ixi[j] * kT;
1763 else /* Other non Trotter temperature NH control -- no chains yet. */
1765 energy += 0.5 * BOLTZ * nd * gmx::square(ivxi[0]) / iQinv[0];
1766 energy += nd * ixi[0] * kT;
1774 /* Returns the energy from the barostat thermostat chain */
1775 static real energyPressureMTTK(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1779 int nh = state->nhchainlength;
1781 for (int i = 0; i < state->nnhpres; i++)
1783 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1784 real reft = std::max<real>(ir->opts.ref_t[0], 0.0); /* using 'System' temperature */
1785 real kT = BOLTZ * reft;
1787 for (int j = 0; j < nh; j++)
1789 double iQinv = MassQ->QPinv[i * nh + j];
1792 energy += 0.5 * gmx::square(state->nhpres_vxi[i * nh + j]) / iQinv;
1793 /* contribution from the thermal variable of the NH chain */
1794 energy += state->nhpres_xi[i * nh + j] * kT;
1799 "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f",
1802 state->nhpres_vxi[i * nh + j],
1803 state->nhpres_xi[i * nh + j]);
1811 /* Returns the energy accumulated by the V-rescale or Berendsen thermostat */
1812 static real energyVrescale(const t_inputrec* ir, const t_state* state)
1815 for (int i = 0; i < ir->opts.ngtc; i++)
1817 energy += state->therm_integral[i];
1823 real NPT_energy(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1827 if (ir->epc != epcNO)
1829 /* Compute the contribution of the pressure to the conserved quantity*/
1831 real vol = det(state->box);
1835 case epcPARRINELLORAHMAN:
1837 /* contribution from the pressure momenta */
1839 calcParrinelloRahmanInvMass(ir, state->box, invMass);
1840 for (int d = 0; d < DIM; d++)
1842 for (int n = 0; n <= d; n++)
1844 if (invMass[d][n] > 0)
1846 energyNPT += 0.5 * gmx::square(state->boxv[d][n]) / (invMass[d][n] * PRESFAC);
1851 /* Contribution from the PV term.
1852 * Not that with non-zero off-diagonal reference pressures,
1853 * i.e. applied shear stresses, there are additional terms.
1854 * We don't support this here, since that requires keeping
1855 * track of unwrapped box diagonal elements. This case is
1856 * excluded in integratorHasConservedEnergyQuantity().
1858 energyNPT += vol * trace(ir->ref_p) / (DIM * PRESFAC);
1862 /* contribution from the pressure momenta */
1863 energyNPT += 0.5 * gmx::square(state->veta) / MassQ->Winv;
1865 /* contribution from the PV term */
1866 energyNPT += vol * trace(ir->ref_p) / (DIM * PRESFAC);
1868 if (ir->epc == epcMTTK)
1870 /* contribution from the MTTK chain */
1871 energyNPT += energyPressureMTTK(ir, state, MassQ);
1875 case epcCRESCALE: energyNPT += state->baros_integral; break;
1879 "Conserved energy quantity for pressure coupling is not handled. A case "
1880 "should be added with either the conserved quantity added or nothing added "
1881 "and an exclusion added to integratorHasConservedEnergyQuantity().");
1889 case etcBERENDSEN: energyNPT += energyVrescale(ir, state); break;
1890 case etcNOSEHOOVER: energyNPT += energyNoseHoover(ir, state, MassQ); break;
1892 case etcANDERSENMASSIVE:
1893 // Not supported, excluded in integratorHasConservedEnergyQuantity()
1898 "Conserved energy quantity for temperature coupling is not handled. A case "
1899 "should be added with either the conserved quantity added or nothing added and "
1900 "an exclusion added to integratorHasConservedEnergyQuantity().");
1907 static real vrescale_sumnoises(real nn, gmx::ThreeFry2x64<>* rng, gmx::NormalDistribution<real>* normalDist)
1910 * Returns the sum of nn independent gaussian noises squared
1911 * (i.e. equivalent to summing the square of the return values
1912 * of nn calls to a normal distribution).
1914 const real ndeg_tol = 0.0001;
1916 gmx::GammaDistribution<real> gammaDist(0.5 * nn, 1.0);
1918 if (nn < 2 + ndeg_tol)
1923 nn_int = gmx::roundToInt(nn);
1925 if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1928 "The v-rescale thermostat was called with a group with #DOF=%f, but for "
1929 "#DOF<3 only integer #DOF are supported",
1934 for (i = 0; i < nn_int; i++)
1936 gauss = (*normalDist)(*rng);
1942 /* Use a gamma distribution for any real nn > 2 */
1943 r = 2.0 * gammaDist(*rng);
1949 real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut, int64_t step, int64_t seed)
1952 * Generates a new value for the kinetic energy,
1953 * according to Bussi et al JCP (2007), Eq. (A7)
1954 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1955 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1956 * ndeg: number of degrees of freedom of the atoms to be thermalized
1957 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1959 real factor, rr, ekin_new;
1960 gmx::ThreeFry2x64<64> rng(seed, gmx::RandomDomain::Thermostat);
1961 gmx::NormalDistribution<real> normalDist;
1965 factor = exp(-1.0 / taut);
1972 rng.restart(step, 0);
1974 rr = normalDist(rng);
1978 * (sigma * (vrescale_sumnoises(ndeg - 1, &rng, &normalDist) + rr * rr) / ndeg - kk)
1979 + 2.0 * rr * std::sqrt(kk * sigma / ndeg * (1.0 - factor) * factor);
1984 void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, real dt, double therm_integral[])
1986 const t_grpopts* opts;
1988 real Ek, Ek_ref1, Ek_ref, Ek_new;
1992 for (i = 0; (i < opts->ngtc); i++)
1996 Ek = trace(ekind->tcstat[i].ekinf);
2000 Ek = trace(ekind->tcstat[i].ekinh);
2003 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
2005 Ek_ref1 = 0.5 * opts->ref_t[i] * BOLTZ;
2006 Ek_ref = Ek_ref1 * opts->nrdf[i];
2008 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i], opts->tau_t[i] / dt, step, ir->ld_seed);
2010 /* Analytically Ek_new>=0, but we check for rounding errors */
2013 ekind->tcstat[i].lambda = 0.0;
2017 ekind->tcstat[i].lambda = std::sqrt(Ek_new / Ek);
2020 therm_integral[i] -= Ek_new - Ek;
2025 "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
2030 ekind->tcstat[i].lambda);
2035 ekind->tcstat[i].lambda = 1.0;
2040 void rescale_velocities(const gmx_ekindata_t* ekind, const t_mdatoms* mdatoms, int start, int end, rvec v[])
2042 const unsigned short* cTC = mdatoms->cTC;
2043 gmx::ArrayRef<const t_grp_tcstat> tcstat = ekind->tcstat;
2045 for (int n = start; n < end; n++)
2052 const real lg = tcstat[gt].lambda;
2053 for (int d = 0; d < DIM; d++)
2060 //! Check whether we're doing simulated annealing
2061 bool doSimulatedAnnealing(const t_inputrec* ir)
2063 for (int i = 0; i < ir->opts.ngtc; i++)
2065 /* set bSimAnn if any group is being annealed */
2066 if (ir->opts.annealing[i] != eannNO)
2074 // TODO If we keep simulated annealing, make a proper module that
2075 // does not rely on changing inputrec.
2076 bool initSimulatedAnnealing(t_inputrec* ir, gmx::Update* upd)
2078 bool doSimAnnealing = doSimulatedAnnealing(ir);
2081 update_annealing_target_temp(ir, ir->init_t, upd);
2083 return doSimAnnealing;
2086 /* set target temperatures if we are annealing */
2087 void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd)
2089 int i, j, n, npoints;
2090 real pert, thist = 0, x;
2092 for (i = 0; i < ir->opts.ngtc; i++)
2094 npoints = ir->opts.anneal_npoints[i];
2095 switch (ir->opts.annealing[i])
2097 case eannNO: continue;
2099 /* calculate time modulo the period */
2100 pert = ir->opts.anneal_time[i][npoints - 1];
2101 n = static_cast<int>(t / pert);
2102 thist = t - n * pert; /* modulo time */
2103 /* Make sure rounding didn't get us outside the interval */
2104 if (std::fabs(thist - pert) < GMX_REAL_EPS * 100)
2109 case eannSINGLE: thist = t; break;
2112 "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",
2117 /* We are doing annealing for this group if we got here,
2118 * and we have the (relative) time as thist.
2119 * calculate target temp */
2121 while ((j < npoints - 1) && (thist > (ir->opts.anneal_time[i][j + 1])))
2125 if (j < npoints - 1)
2127 /* Found our position between points j and j+1.
2128 * Interpolate: x is the amount from j+1, (1-x) from point j
2129 * First treat possible jumps in temperature as a special case.
2131 if ((ir->opts.anneal_time[i][j + 1] - ir->opts.anneal_time[i][j]) < GMX_REAL_EPS * 100)
2133 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][j + 1];
2137 x = ((thist - ir->opts.anneal_time[i][j])
2138 / (ir->opts.anneal_time[i][j + 1] - ir->opts.anneal_time[i][j]));
2140 x * ir->opts.anneal_temp[i][j + 1] + (1 - x) * ir->opts.anneal_temp[i][j];
2145 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][npoints - 1];
2149 upd->update_temperature_constants(*ir);
2152 void pleaseCiteCouplingAlgorithms(FILE* fplog, const t_inputrec& ir)
2154 if (EI_DYNAMICS(ir.eI))
2156 if (ir.etc == etcBERENDSEN)
2158 please_cite(fplog, "Berendsen84a");
2160 if (ir.etc == etcVRESCALE)
2162 please_cite(fplog, "Bussi2007a");
2164 if (ir.epc == epcCRESCALE)
2166 please_cite(fplog, "Bernetti2020");
2168 // TODO this is actually an integrator, not a coupling algorithm
2171 please_cite(fplog, "Goga2012");