2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
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, 0.2967324292201065, -0.186929716880426,
91 0.2967324292201065, 0.2967324292201065 };
93 static const double* sy_const[] = { nullptr, sy_const_1, nullptr, sy_const_3, nullptr, sy_const_5 };
96 void update_tcouple(int64_t step,
97 const t_inputrec* inputrec,
99 gmx_ekindata_t* ekind,
100 const t_extmass* MassQ,
104 // This condition was explicitly checked in previous version, but should have never been satisfied
105 GMX_ASSERT(!(EI_VV(inputrec->eI)
106 && (inputrecNvtTrotter(inputrec) || inputrecNptTrotter(inputrec)
107 || inputrecNphTrotter(inputrec))),
108 "Temperature coupling was requested with velocity verlet and trotter");
110 bool doTemperatureCoupling = false;
112 // For VV temperature coupling parameters are updated on the current
113 // step, for the others - one step before.
114 if (inputrec->etc == etcNO)
116 doTemperatureCoupling = false;
118 else if (EI_VV(inputrec->eI))
120 doTemperatureCoupling = do_per_step(step, inputrec->nsttcouple);
124 doTemperatureCoupling = do_per_step(step + inputrec->nsttcouple - 1, inputrec->nsttcouple);
127 if (doTemperatureCoupling)
129 real dttc = inputrec->nsttcouple * inputrec->delta_t;
131 // TODO: berendsen_tcoupl(...), nosehoover_tcoupl(...) and vrescale_tcoupl(...) update
132 // temperature coupling parameters, which should be reflected in the name of these
134 switch (inputrec->etc)
138 berendsen_tcoupl(inputrec, ekind, dttc, state->therm_integral);
141 nosehoover_tcoupl(&(inputrec->opts), ekind, dttc, state->nosehoover_xi.data(),
142 state->nosehoover_vxi.data(), MassQ);
145 vrescale_tcoupl(inputrec, step, ekind, dttc, state->therm_integral.data());
148 /* rescale in place here */
149 if (EI_VV(inputrec->eI))
151 rescale_velocities(ekind, md, 0, md->homenr, state->v.rvec_array());
156 // Set the T scaling lambda to 1 to have no scaling
157 // TODO: Do we have to do it on every non-t-couple step?
158 for (int i = 0; (i < inputrec->opts.ngtc); i++)
160 ekind->tcstat[i].lambda = 1.0;
165 void update_pcouple_before_coordinates(FILE* fplog,
167 const t_inputrec* inputrec,
169 matrix parrinellorahmanMu,
173 /* Berendsen P-coupling is completely handled after the coordinate update.
174 * Trotter P-coupling is handled by separate calls to trotter_update().
176 if (inputrec->epc == epcPARRINELLORAHMAN
177 && do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
179 real dtpc = inputrec->nstpcouple * inputrec->delta_t;
181 parrinellorahman_pcoupl(fplog, step, inputrec, dtpc, state->pres_prev, state->box,
182 state->box_rel, state->boxv, M, parrinellorahmanMu, bInitStep);
186 void update_pcouple_after_coordinates(FILE* fplog,
188 const t_inputrec* inputrec,
190 const matrix pressure,
191 const matrix forceVirial,
192 const matrix constraintVirial,
193 matrix pressureCouplingMu,
196 gmx::BoxDeformation* boxDeformation,
197 const bool scaleCoordinates)
200 int homenr = md->homenr;
202 /* Cast to real for faster code, no loss in precision (see comment above) */
203 real dt = inputrec->delta_t;
206 /* now update boxes */
207 switch (inputrec->epc)
211 if (do_per_step(step, inputrec->nstpcouple))
213 real dtpc = inputrec->nstpcouple * dt;
214 berendsen_pcoupl(fplog, step, inputrec, dtpc, pressure, state->box, forceVirial,
215 constraintVirial, pressureCouplingMu, &state->baros_integral);
216 berendsen_pscale(inputrec, pressureCouplingMu, state->box, state->box_rel, start,
217 homenr, state->x.rvec_array(), md->cFREEZE, nrnb, scaleCoordinates);
220 case (epcPARRINELLORAHMAN):
221 if (do_per_step(step + inputrec->nstpcouple - 1, inputrec->nstpcouple))
223 /* The box velocities were updated in do_pr_pcoupl,
224 * but we dont change the box vectors until we get here
225 * since we need to be able to shift/unshift above.
227 real dtpc = inputrec->nstpcouple * dt;
228 for (int i = 0; i < DIM; i++)
230 for (int m = 0; m <= i; m++)
232 state->box[i][m] += dtpc * state->boxv[i][m];
235 preserve_box_shape(inputrec, state->box_rel, state->box);
237 /* Scale the coordinates */
238 if (scaleCoordinates)
240 auto x = state->x.rvec_array();
241 for (int n = start; n < start + homenr; n++)
243 tmvmul_ur0(pressureCouplingMu, x[n], x[n]);
249 switch (inputrec->epct)
251 case (epctISOTROPIC):
252 /* DIM * eta = ln V. so DIM*eta_new = DIM*eta_old + DIM*dt*veta =>
253 ln V_new = ln V_old + 3*dt*veta => V_new = V_old*exp(3*dt*veta) =>
254 Side length scales as exp(veta*dt) */
256 msmul(state->box, std::exp(state->veta * dt), state->box);
258 /* Relate veta to boxv. veta = d(eta)/dT = (1/DIM)*1/V dV/dT.
259 o If we assume isotropic scaling, and box length scaling
260 factor L, then V = L^DIM (det(M)). So dV/dt = DIM
261 L^(DIM-1) dL/dt det(M), and veta = (1/L) dL/dt. The
262 determinant of B is L^DIM det(M), and the determinant
263 of dB/dt is (dL/dT)^DIM det (M). veta will be
264 (det(dB/dT)/det(B))^(1/3). Then since M =
265 B_new*(vol_new)^(1/3), dB/dT_new = (veta_new)*B(new). */
267 msmul(state->box, state->veta, state->boxv);
277 auto localX = makeArrayRef(state->x).subArray(start, homenr);
278 boxDeformation->apply(localX, state->box, step);
282 extern gmx_bool update_randomize_velocities(const t_inputrec* ir,
286 gmx::ArrayRef<gmx::RVec> v,
287 const gmx::Update* upd,
288 const gmx::Constraints* constr)
291 real rate = (ir->delta_t) / ir->opts.tau_t[0];
293 if (ir->etc == etcANDERSEN && constr != nullptr)
295 /* Currently, Andersen thermostat does not support constrained
296 systems. Functionality exists in the andersen_tcoupl
297 function in GROMACS 4.5.7 to allow this combination. That
298 code could be ported to the current random-number
299 generation approach, but has not yet been done because of
300 lack of time and resources. */
302 "Normal Andersen is currently not supported with constraints, use massive "
306 /* proceed with andersen if 1) it's fixed probability per
307 particle andersen or 2) it's massive andersen and it's tau_t/dt */
308 if ((ir->etc == etcANDERSEN) || do_per_step(step, gmx::roundToInt(1.0 / rate)))
310 andersen_tcoupl(ir, step, cr, md, v, rate, upd->getAndersenRandomizeGroup(),
311 upd->getBoltzmanFactor());
318 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
322 {0.828981543588751,-0.657963087177502,0.828981543588751},
324 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
327 /* these integration routines are only referenced inside this file */
328 static void NHC_trotter(const t_grpopts* opts,
330 const gmx_ekindata_t* ekind,
336 const t_extmass* MassQ,
337 gmx_bool bEkinAveVel)
340 /* general routine for both barostat and thermostat nose hoover chains */
343 double Ekin, Efac, reft, kT, nd;
348 int mstepsi, mstepsj;
349 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
350 int nh = opts->nhchainlength;
353 mstepsi = mstepsj = ns;
355 /* if scalefac is NULL, we are doing the NHC of the barostat */
358 if (scalefac == nullptr)
363 for (i = 0; i < nvar; i++)
366 /* make it easier to iterate by selecting
367 out the sub-array that corresponds to this T group */
371 gmx::ArrayRef<const double> iQinv;
374 iQinv = gmx::arrayRefFromArray(&MassQ->QPinv[i * nh], nh);
375 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
376 reft = std::max<real>(0, opts->ref_t[0]);
377 Ekin = gmx::square(*veta) / MassQ->Winv;
381 iQinv = gmx::arrayRefFromArray(&MassQ->Qinv[i * nh], nh);
382 const t_grp_tcstat* tcstat = &ekind->tcstat[i];
384 reft = std::max<real>(0, opts->ref_t[i]);
387 Ekin = 2 * trace(tcstat->ekinf) * tcstat->ekinscalef_nhc;
391 Ekin = 2 * trace(tcstat->ekinh) * tcstat->ekinscaleh_nhc;
396 for (mi = 0; mi < mstepsi; mi++)
398 for (mj = 0; mj < mstepsj; mj++)
400 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
401 dt = sy_const[ns][mj] * dtfull / mstepsi;
403 /* compute the thermal forces */
404 GQ[0] = iQinv[0] * (Ekin - nd * kT);
406 for (j = 0; j < nh - 1; j++)
408 if (iQinv[j + 1] > 0)
410 /* we actually don't need to update here if we save the
411 state of the GQ, but it's easier to just recompute*/
412 GQ[j + 1] = iQinv[j + 1] * ((gmx::square(ivxi[j]) / iQinv[j]) - kT);
420 ivxi[nh - 1] += 0.25 * dt * GQ[nh - 1];
421 for (j = nh - 1; j > 0; j--)
423 Efac = exp(-0.125 * dt * ivxi[j]);
424 ivxi[j - 1] = Efac * (ivxi[j - 1] * Efac + 0.25 * dt * GQ[j - 1]);
427 Efac = exp(-0.5 * dt * ivxi[0]);
436 Ekin *= (Efac * Efac);
438 /* Issue - if the KE is an average of the last and the current temperatures, then we
439 might not be able to scale the kinetic energy directly with this factor. Might
440 take more bookkeeping -- have to think about this a bit more . . . */
442 GQ[0] = iQinv[0] * (Ekin - nd * kT);
444 /* update thermostat positions */
445 for (j = 0; j < nh; j++)
447 ixi[j] += 0.5 * dt * ivxi[j];
450 for (j = 0; j < nh - 1; j++)
452 Efac = exp(-0.125 * dt * ivxi[j + 1]);
453 ivxi[j] = Efac * (ivxi[j] * Efac + 0.25 * dt * GQ[j]);
454 if (iQinv[j + 1] > 0)
456 GQ[j + 1] = iQinv[j + 1] * ((gmx::square(ivxi[j]) / iQinv[j]) - kT);
463 ivxi[nh - 1] += 0.25 * dt * GQ[nh - 1];
470 static void boxv_trotter(const t_inputrec* ir,
474 const gmx_ekindata_t* ekind,
477 const t_extmass* MassQ)
484 tensor ekinmod, localpres;
486 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
487 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
490 if (ir->epct == epctSEMIISOTROPIC)
499 /* eta is in pure units. veta is in units of ps^-1. GW is in
500 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
501 taken to use only RATIOS of eta in updating the volume. */
503 /* we take the partial pressure tensors, modify the
504 kinetic energy tensor, and recovert to pressure */
506 if (ir->opts.nrdf[0] == 0)
508 gmx_fatal(FARGS, "Barostat is coupled to a T-group with no degrees of freedom\n");
510 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
511 alpha = 1.0 + DIM / (static_cast<double>(ir->opts.nrdf[0]));
512 alpha *= ekind->tcstat[0].ekinscalef_nhc;
513 msmul(ekind->ekin, alpha, ekinmod);
514 /* for now, we use Elr = 0, because if you want to get it right, you
515 really should be using PME. Maybe print a warning? */
517 pscal = calc_pres(ir->pbcType, nwall, box, ekinmod, vir, localpres) + pcorr;
520 GW = (vol * (MassQ->Winv / PRESFAC)) * (DIM * pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
522 *veta += 0.5 * dt * GW;
526 * This file implements temperature and pressure coupling algorithms:
527 * For now only the Weak coupling and the modified weak coupling.
529 * Furthermore computation of pressure and temperature is done here
533 real calc_pres(PbcType pbcType, int nwall, const matrix box, const tensor ekin, const tensor vir, tensor pres)
538 if (pbcType == PbcType::No || (pbcType == PbcType::XY && nwall != 2))
544 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
545 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
549 fac = PRESFAC * 2.0 / det(box);
550 for (n = 0; (n < DIM); n++)
552 for (m = 0; (m < DIM); m++)
554 pres[n][m] = (ekin[n][m] - vir[n][m]) * fac;
560 pr_rvecs(debug, 0, "PC: pres", pres, DIM);
561 pr_rvecs(debug, 0, "PC: ekin", ekin, DIM);
562 pr_rvecs(debug, 0, "PC: vir ", vir, DIM);
563 pr_rvecs(debug, 0, "PC: box ", box, DIM);
566 return trace(pres) / DIM;
569 real calc_temp(real ekin, real nrdf)
573 return (2.0 * ekin) / (nrdf * BOLTZ);
581 /*! \brief Sets 1/mass for Parrinello-Rahman in wInv; NOTE: PRESFAC is not included, so not in GROMACS units! */
582 static void calcParrinelloRahmanInvMass(const t_inputrec* ir, const matrix box, tensor wInv)
586 /* TODO: See if we can make the mass independent of the box size */
587 maxBoxLength = std::max(box[XX][XX], box[YY][YY]);
588 maxBoxLength = std::max(maxBoxLength, box[ZZ][ZZ]);
590 for (int d = 0; d < DIM; d++)
592 for (int n = 0; n < DIM; n++)
594 wInv[d][n] = (4 * M_PI * M_PI * ir->compress[d][n]) / (3 * ir->tau_p * ir->tau_p * maxBoxLength);
599 void parrinellorahman_pcoupl(FILE* fplog,
601 const t_inputrec* ir,
611 /* This doesn't do any coordinate updating. It just
612 * integrates the box vector equations from the calculated
613 * acceleration due to pressure difference. We also compute
614 * the tensor M which is used in update to couple the particle
615 * coordinates to the box vectors.
617 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
620 * M_nk = (h') * (h' * h + h' h) * h
622 * with the dots denoting time derivatives and h is the transformation from
623 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
624 * This also goes for the pressure and M tensors - they are transposed relative
625 * to ours. Our equation thus becomes:
628 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
630 * where b is the gromacs box matrix.
631 * Our box accelerations are given by
633 * b = vol/W inv(box') * (P-ref_P) (=h')
636 real vol = box[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
637 real atot, arel, change;
638 tensor invbox, pdiff, t1, t2;
640 gmx::invertBoxMatrix(box, invbox);
644 /* Note that PRESFAC does not occur here.
645 * The pressure and compressibility always occur as a product,
646 * therefore the pressure unit drops out.
649 calcParrinelloRahmanInvMass(ir, box, winv);
651 m_sub(pres, ir->ref_p, pdiff);
653 if (ir->epct == epctSURFACETENSION)
655 /* Unlike Berendsen coupling it might not be trivial to include a z
656 * pressure correction here? On the other hand we don't scale the
657 * box momentarily, but change accelerations, so it might not be crucial.
659 real xy_pressure = 0.5 * (pres[XX][XX] + pres[YY][YY]);
660 for (int d = 0; d < ZZ; d++)
662 pdiff[d][d] = (xy_pressure - (pres[ZZ][ZZ] - ir->ref_p[d][d] / box[d][d]));
666 tmmul(invbox, pdiff, t1);
667 /* Move the off-diagonal elements of the 'force' to one side to ensure
668 * that we obey the box constraints.
670 for (int d = 0; d < DIM; d++)
672 for (int n = 0; n < d; n++)
674 t1[d][n] += t1[n][d];
681 case epctANISOTROPIC:
682 for (int d = 0; d < DIM; d++)
684 for (int n = 0; n <= d; n++)
686 t1[d][n] *= winv[d][n] * vol;
691 /* calculate total volume acceleration */
692 atot = box[XX][XX] * box[YY][YY] * t1[ZZ][ZZ] + box[XX][XX] * t1[YY][YY] * box[ZZ][ZZ]
693 + t1[XX][XX] * box[YY][YY] * box[ZZ][ZZ];
694 arel = atot / (3 * vol);
695 /* set all RELATIVE box accelerations equal, and maintain total V
697 for (int d = 0; d < DIM; d++)
699 for (int n = 0; n <= d; n++)
701 t1[d][n] = winv[0][0] * vol * arel * box[d][n];
705 case epctSEMIISOTROPIC:
706 case epctSURFACETENSION:
707 /* Note the correction to pdiff above for surftens. coupling */
709 /* calculate total XY volume acceleration */
710 atot = box[XX][XX] * t1[YY][YY] + t1[XX][XX] * box[YY][YY];
711 arel = atot / (2 * box[XX][XX] * box[YY][YY]);
712 /* set RELATIVE XY box accelerations equal, and maintain total V
713 * change speed. Dont change the third box vector accelerations */
714 for (int d = 0; d < ZZ; d++)
716 for (int n = 0; n <= d; n++)
718 t1[d][n] = winv[d][n] * vol * arel * box[d][n];
721 for (int n = 0; n < DIM; n++)
723 t1[ZZ][n] *= winv[ZZ][n] * vol;
728 "Parrinello-Rahman pressure coupling type %s "
729 "not supported yet\n",
730 EPCOUPLTYPETYPE(ir->epct));
734 for (int d = 0; d < DIM; d++)
736 for (int n = 0; n <= d; n++)
738 boxv[d][n] += dt * t1[d][n];
740 /* We do NOT update the box vectors themselves here, since
741 * we need them for shifting later. It is instead done last
742 * in the update() routine.
745 /* Calculate the change relative to diagonal elements-
746 since it's perfectly ok for the off-diagonal ones to
747 be zero it doesn't make sense to check the change relative
751 change = std::fabs(dt * boxv[d][n] / box[d][d]);
753 if (change > maxchange)
760 if (maxchange > 0.01 && fplog)
764 "\nStep %s Warning: Pressure scaling more than 1%%. "
765 "This may mean your system\n is not yet equilibrated. "
766 "Use of Parrinello-Rahman pressure coupling during\n"
767 "equilibration can lead to simulation instability, "
768 "and is discouraged.\n",
769 gmx_step_str(step, buf));
773 preserve_box_shape(ir, box_rel, boxv);
775 mtmul(boxv, box, t1); /* t1=boxv * b' */
776 mmul(invbox, t1, t2);
777 mtmul(t2, invbox, M);
779 /* Determine the scaling matrix mu for the coordinates */
780 for (int d = 0; d < DIM; d++)
782 for (int n = 0; n <= d; n++)
784 t1[d][n] = box[d][n] + dt * boxv[d][n];
787 preserve_box_shape(ir, box_rel, t1);
788 /* t1 is the box at t+dt, determine mu as the relative change */
789 mmul_ur0(invbox, t1, mu);
792 void berendsen_pcoupl(FILE* fplog,
794 const t_inputrec* ir,
798 const matrix force_vir,
799 const matrix constraint_vir,
801 double* baros_integral)
804 real scalar_pressure, xy_pressure, p_corr_z;
808 * Calculate the scaling matrix mu
812 for (d = 0; d < DIM; d++)
814 scalar_pressure += pres[d][d] / DIM;
817 xy_pressure += pres[d][d] / (DIM - 1);
820 /* Pressure is now in bar, everywhere. */
821 #define factor(d, m) (ir->compress[d][m] * dt / ir->tau_p)
823 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
824 * necessary for triclinic scaling
830 for (d = 0; d < DIM; d++)
832 mu[d][d] = 1.0 - factor(d, d) * (ir->ref_p[d][d] - scalar_pressure) / DIM;
835 case epctSEMIISOTROPIC:
836 for (d = 0; d < ZZ; d++)
838 mu[d][d] = 1.0 - factor(d, d) * (ir->ref_p[d][d] - xy_pressure) / DIM;
840 mu[ZZ][ZZ] = 1.0 - factor(ZZ, ZZ) * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]) / DIM;
842 case epctANISOTROPIC:
843 for (d = 0; d < DIM; d++)
845 for (n = 0; n < DIM; n++)
847 mu[d][n] = (d == n ? 1.0 : 0.0) - factor(d, n) * (ir->ref_p[d][n] - pres[d][n]) / DIM;
851 case epctSURFACETENSION:
852 /* ir->ref_p[0/1] is the reference surface-tension times *
853 * the number of surfaces */
854 if (ir->compress[ZZ][ZZ] != 0.0F)
856 p_corr_z = dt / ir->tau_p * (ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
860 /* when the compressibity is zero, set the pressure correction *
861 * in the z-direction to zero to get the correct surface tension */
864 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ] * p_corr_z;
865 for (d = 0; d < DIM - 1; d++)
869 * (ir->ref_p[d][d] / (mu[ZZ][ZZ] * box[ZZ][ZZ])
870 - (pres[ZZ][ZZ] + p_corr_z - xy_pressure))
875 gmx_fatal(FARGS, "Berendsen pressure coupling type %s not supported yet\n",
876 EPCOUPLTYPETYPE(ir->epct));
878 /* To fullfill the orientation restrictions on triclinic boxes
879 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
880 * the other elements of mu to first order.
882 mu[YY][XX] += mu[XX][YY];
883 mu[ZZ][XX] += mu[XX][ZZ];
884 mu[ZZ][YY] += mu[YY][ZZ];
889 /* Keep track of the work the barostat applies on the system.
890 * Without constraints force_vir tells us how Epot changes when scaling.
891 * With constraints constraint_vir gives us the constraint contribution
892 * to both Epot and Ekin. Although we are not scaling velocities, scaling
893 * the coordinates leads to scaling of distances involved in constraints.
894 * This in turn changes the angular momentum (even if the constrained
895 * distances are corrected at the next step). The kinetic component
896 * of the constraint virial captures the angular momentum change.
898 for (int d = 0; d < DIM; d++)
900 for (int n = 0; n <= d; n++)
903 2 * (mu[d][n] - (n == d ? 1 : 0)) * (force_vir[d][n] + constraint_vir[d][n]);
909 pr_rvecs(debug, 0, "PC: pres ", pres, 3);
910 pr_rvecs(debug, 0, "PC: mu ", mu, 3);
913 if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 || mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01
914 || mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
918 "\nStep %s Warning: pressure scaling more than 1%%, "
920 gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
923 fprintf(fplog, "%s", buf);
925 fprintf(stderr, "%s", buf);
929 void berendsen_pscale(const t_inputrec* ir,
936 const unsigned short cFREEZE[],
938 const bool scaleCoordinates)
940 ivec* nFreeze = ir->opts.nFreeze;
942 int nthreads gmx_unused;
944 #ifndef __clang_analyzer__
945 nthreads = gmx_omp_nthreads_get(emntUpdate);
948 /* Scale the positions */
949 if (scaleCoordinates)
951 #pragma omp parallel for num_threads(nthreads) schedule(static)
952 for (int n = start; n < start + nr_atoms; n++)
954 // Trivial OpenMP region that does not throw
957 if (cFREEZE == nullptr)
968 x[n][XX] = mu[XX][XX] * x[n][XX] + mu[YY][XX] * x[n][YY] + mu[ZZ][XX] * x[n][ZZ];
972 x[n][YY] = mu[YY][YY] * x[n][YY] + mu[ZZ][YY] * x[n][ZZ];
976 x[n][ZZ] = mu[ZZ][ZZ] * x[n][ZZ];
980 /* compute final boxlengths */
981 for (d = 0; d < DIM; d++)
983 box[d][XX] = mu[XX][XX] * box[d][XX] + mu[YY][XX] * box[d][YY] + mu[ZZ][XX] * box[d][ZZ];
984 box[d][YY] = mu[YY][YY] * box[d][YY] + mu[ZZ][YY] * box[d][ZZ];
985 box[d][ZZ] = mu[ZZ][ZZ] * box[d][ZZ];
988 preserve_box_shape(ir, box_rel, box);
990 /* (un)shifting should NOT be done after this,
991 * since the box vectors might have changed
993 inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
996 void berendsen_tcoupl(const t_inputrec* ir, gmx_ekindata_t* ekind, real dt, std::vector<double>& therm_integral)
998 const t_grpopts* opts = &ir->opts;
1000 for (int i = 0; (i < opts->ngtc); i++)
1006 Ek = trace(ekind->tcstat[i].ekinf);
1007 T = ekind->tcstat[i].T;
1011 Ek = trace(ekind->tcstat[i].ekinh);
1012 T = ekind->tcstat[i].Th;
1015 if ((opts->tau_t[i] > 0) && (T > 0.0))
1017 real reft = std::max<real>(0, opts->ref_t[i]);
1018 real lll = std::sqrt(1.0 + (dt / opts->tau_t[i]) * (reft / T - 1.0));
1019 ekind->tcstat[i].lambda = std::max<real>(std::min<real>(lll, 1.25), 0.8);
1023 ekind->tcstat[i].lambda = 1.0;
1026 /* Keep track of the amount of energy we are adding to the system */
1027 therm_integral[i] -= (gmx::square(ekind->tcstat[i].lambda) - 1) * Ek;
1031 fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n", i, T, ekind->tcstat[i].lambda);
1036 void andersen_tcoupl(const t_inputrec* ir,
1038 const t_commrec* cr,
1039 const t_mdatoms* md,
1040 gmx::ArrayRef<gmx::RVec> v,
1042 const std::vector<bool>& randomize,
1043 gmx::ArrayRef<const real> boltzfac)
1045 const int* gatindex = (DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
1048 gmx::ThreeFry2x64<0> rng(ir->andersen_seed, gmx::RandomDomain::Thermostat);
1049 gmx::UniformRealDistribution<real> uniformDist;
1050 gmx::TabulatedNormalDistribution<real, 14> normalDist;
1052 /* randomize the velocities of the selected particles */
1054 for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */
1056 int ng = gatindex ? gatindex[i] : i;
1057 gmx_bool bRandomize;
1059 rng.restart(step, ng);
1063 gc = md->cTC[i]; /* assign the atom to a temperature group if there are more than one */
1067 if (ir->etc == etcANDERSENMASSIVE)
1069 /* Randomize particle always */
1074 /* Randomize particle probabilistically */
1075 uniformDist.reset();
1076 bRandomize = uniformDist(rng) < rate;
1083 scal = std::sqrt(boltzfac[gc] * md->invmass[i]);
1087 for (d = 0; d < DIM; d++)
1089 v[i][d] = scal * normalDist(rng);
1097 void nosehoover_tcoupl(const t_grpopts* opts,
1098 const gmx_ekindata_t* ekind,
1102 const t_extmass* MassQ)
1107 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
1109 for (i = 0; (i < opts->ngtc); i++)
1111 reft = std::max<real>(0, opts->ref_t[i]);
1113 vxi[i] += dt * MassQ->Qinv[i] * (ekind->tcstat[i].Th - reft);
1114 xi[i] += dt * (oldvxi + vxi[i]) * 0.5;
1118 void trotter_update(const t_inputrec* ir,
1120 gmx_ekindata_t* ekind,
1121 const gmx_enerdata_t* enerd,
1124 const t_mdatoms* md,
1125 const t_extmass* MassQ,
1126 gmx::ArrayRef<std::vector<int>> trotter_seqlist,
1130 int n, i, d, ngtc, gc = 0, t;
1131 t_grp_tcstat* tcstat;
1132 const t_grpopts* opts;
1135 double * scalefac, dtc;
1136 rvec sumv = { 0, 0, 0 };
1139 if (trotter_seqno <= ettTSEQ2)
1141 step_eff = step - 1; /* the velocity verlet calls are actually out of order -- the first
1142 half step is actually the last half step from the previous step.
1143 Thus the first half step actually corresponds to the n-1 step*/
1150 bCouple = (ir->nsttcouple == 1 || do_per_step(step_eff + ir->nsttcouple, ir->nsttcouple));
1152 const gmx::ArrayRef<const int> trotter_seq = trotter_seqlist[trotter_seqno];
1154 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
1158 dtc = ir->nsttcouple * ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
1159 opts = &(ir->opts); /* just for ease of referencing */
1162 snew(scalefac, opts->ngtc);
1163 for (i = 0; i < ngtc; i++)
1167 /* execute the series of trotter updates specified in the trotterpart array */
1169 for (i = 0; i < NTROTTERPARTS; i++)
1171 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
1172 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2)
1173 || (trotter_seq[i] == etrtNHC2))
1182 auto v = makeArrayRef(state->v);
1183 switch (trotter_seq[i])
1187 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir,
1188 enerd->term[F_PDISPCORR], MassQ);
1192 NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi.data(),
1193 state->nhpres_vxi.data(), nullptr, &(state->veta), MassQ, FALSE);
1197 NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi.data(),
1198 state->nosehoover_vxi.data(), scalefac, nullptr, MassQ, (ir->eI == eiVV));
1199 /* need to rescale the kinetic energies and velocities here. Could
1200 scale the velocities later, but we need them scaled in order to
1201 produce the correct outputs, so we'll scale them here. */
1203 for (t = 0; t < ngtc; t++)
1205 tcstat = &ekind->tcstat[t];
1206 tcstat->vscale_nhc = scalefac[t];
1207 tcstat->ekinscaleh_nhc *= (scalefac[t] * scalefac[t]);
1208 tcstat->ekinscalef_nhc *= (scalefac[t] * scalefac[t]);
1210 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
1211 /* but do we actually need the total? */
1213 /* modify the velocities as well */
1214 for (n = 0; n < md->homenr; n++)
1216 if (md->cTC) /* does this conditional need to be here? is this always true?*/
1220 for (d = 0; d < DIM; d++)
1222 v[n][d] *= scalefac[gc];
1227 for (d = 0; d < DIM; d++)
1229 sumv[d] += (v[n][d]) / md->invmass[n];
1237 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
1242 extern void init_npt_masses(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool bInit)
1244 int n, i, j, d, ngtc, nh;
1245 const t_grpopts* opts;
1246 real reft, kT, ndj, nd;
1248 opts = &(ir->opts); /* just for ease of referencing */
1249 ngtc = ir->opts.ngtc;
1250 nh = state->nhchainlength;
1256 MassQ->Qinv.resize(ngtc);
1258 for (i = 0; (i < ngtc); i++)
1260 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1262 MassQ->Qinv[i] = 1.0 / (gmx::square(opts->tau_t[i] / M_2PI) * opts->ref_t[i]);
1266 MassQ->Qinv[i] = 0.0;
1270 else if (EI_VV(ir->eI))
1272 /* Set pressure variables */
1276 if (state->vol0 == 0)
1278 state->vol0 = det(state->box);
1279 /* because we start by defining a fixed
1280 compressibility, we need the volume at this
1281 compressibility to solve the problem. */
1285 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1286 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1287 MassQ->Winv = (PRESFAC * trace(ir->compress) * BOLTZ * opts->ref_t[0])
1288 / (DIM * state->vol0 * gmx::square(ir->tau_p / M_2PI));
1289 /* An alternate mass definition, from Tuckerman et al. */
1290 /* MassQ->Winv = 1.0/(gmx::square(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1291 for (d = 0; d < DIM; d++)
1293 for (n = 0; n < DIM; n++)
1295 MassQ->Winvm[d][n] =
1296 PRESFAC * ir->compress[d][n] / (state->vol0 * gmx::square(ir->tau_p / M_2PI));
1297 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1298 before using MTTK for anisotropic states.*/
1301 /* Allocate space for thermostat variables */
1304 MassQ->Qinv.resize(ngtc * nh);
1307 /* now, set temperature variables */
1308 for (i = 0; i < ngtc; i++)
1310 if (opts->tau_t[i] > 0 && opts->ref_t[i] > 0 && opts->nrdf[i] > 0)
1312 reft = std::max<real>(0, opts->ref_t[i]);
1315 for (j = 0; j < nh; j++)
1325 MassQ->Qinv[i * nh + j] = 1.0 / (gmx::square(opts->tau_t[i] / M_2PI) * ndj * kT);
1330 for (j = 0; j < nh; j++)
1332 MassQ->Qinv[i * nh + j] = 0.0;
1339 std::array<std::vector<int>, ettTSEQMAX>
1340 init_npt_vars(const t_inputrec* ir, t_state* state, t_extmass* MassQ, gmx_bool bTrotter)
1342 int i, j, nnhpres, nh;
1343 const t_grpopts* opts;
1344 real bmass, qmass, reft, kT;
1346 opts = &(ir->opts); /* just for ease of referencing */
1347 nnhpres = state->nnhpres;
1348 nh = state->nhchainlength;
1350 if (EI_VV(ir->eI) && (ir->epc == epcMTTK) && (ir->etc != etcNOSEHOOVER))
1353 "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
1356 init_npt_masses(ir, state, MassQ, TRUE);
1358 /* first, initialize clear all the trotter calls */
1359 std::array<std::vector<int>, ettTSEQMAX> trotter_seq;
1360 for (i = 0; i < ettTSEQMAX; i++)
1362 trotter_seq[i].resize(NTROTTERPARTS, etrtNONE);
1363 trotter_seq[i][0] = etrtSKIPALL;
1368 /* no trotter calls, so we never use the values in the array.
1369 * We access them (so we need to define them, but ignore
1375 /* compute the kinetic energy by using the half step velocities or
1376 * the kinetic energies, depending on the order of the trotter calls */
1380 if (inputrecNptTrotter(ir))
1382 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1383 We start with the initial one. */
1384 /* first, a round that estimates veta. */
1385 trotter_seq[0][0] = etrtBAROV;
1387 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1389 /* The first half trotter update */
1390 trotter_seq[2][0] = etrtBAROV;
1391 trotter_seq[2][1] = etrtNHC;
1392 trotter_seq[2][2] = etrtBARONHC;
1394 /* The second half trotter update */
1395 trotter_seq[3][0] = etrtBARONHC;
1396 trotter_seq[3][1] = etrtNHC;
1397 trotter_seq[3][2] = etrtBAROV;
1399 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1401 else if (inputrecNvtTrotter(ir))
1403 /* This is the easy version - there are only two calls, both the same.
1404 Otherwise, even easier -- no calls */
1405 trotter_seq[2][0] = etrtNHC;
1406 trotter_seq[3][0] = etrtNHC;
1408 else if (inputrecNphTrotter(ir))
1410 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1411 We start with the initial one. */
1412 /* first, a round that estimates veta. */
1413 trotter_seq[0][0] = etrtBAROV;
1415 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1417 /* The first half trotter update */
1418 trotter_seq[2][0] = etrtBAROV;
1419 trotter_seq[2][1] = etrtBARONHC;
1421 /* The second half trotter update */
1422 trotter_seq[3][0] = etrtBARONHC;
1423 trotter_seq[3][1] = etrtBAROV;
1425 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1428 else if (ir->eI == eiVVAK)
1430 if (inputrecNptTrotter(ir))
1432 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1433 We start with the initial one. */
1434 /* first, a round that estimates veta. */
1435 trotter_seq[0][0] = etrtBAROV;
1437 /* The first half trotter update, part 1 -- double update, because it commutes */
1438 trotter_seq[1][0] = etrtNHC;
1440 /* The first half trotter update, part 2 */
1441 trotter_seq[2][0] = etrtBAROV;
1442 trotter_seq[2][1] = etrtBARONHC;
1444 /* The second half trotter update, part 1 */
1445 trotter_seq[3][0] = etrtBARONHC;
1446 trotter_seq[3][1] = etrtBAROV;
1448 /* The second half trotter update */
1449 trotter_seq[4][0] = etrtNHC;
1451 else if (inputrecNvtTrotter(ir))
1453 /* This is the easy version - there is only one call, both the same.
1454 Otherwise, even easier -- no calls */
1455 trotter_seq[1][0] = etrtNHC;
1456 trotter_seq[4][0] = etrtNHC;
1458 else if (inputrecNphTrotter(ir))
1460 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1461 We start with the initial one. */
1462 /* first, a round that estimates veta. */
1463 trotter_seq[0][0] = etrtBAROV;
1465 /* The first half trotter update, part 1 -- leave zero */
1466 trotter_seq[1][0] = etrtNHC;
1468 /* The first half trotter update, part 2 */
1469 trotter_seq[2][0] = etrtBAROV;
1470 trotter_seq[2][1] = etrtBARONHC;
1472 /* The second half trotter update, part 1 */
1473 trotter_seq[3][0] = etrtBARONHC;
1474 trotter_seq[3][1] = etrtBAROV;
1476 /* The second half trotter update -- blank for now */
1483 default: bmass = DIM * DIM; /* recommended mass parameters for isotropic barostat */
1486 MassQ->QPinv.resize(nnhpres * opts->nhchainlength);
1488 /* barostat temperature */
1489 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1491 reft = std::max<real>(0, opts->ref_t[0]);
1493 for (i = 0; i < nnhpres; i++)
1495 for (j = 0; j < nh; j++)
1505 MassQ->QPinv[i * opts->nhchainlength + j] =
1506 1.0 / (gmx::square(opts->tau_t[0] / M_2PI) * qmass * kT);
1512 for (i = 0; i < nnhpres; i++)
1514 for (j = 0; j < nh; j++)
1516 MassQ->QPinv[i * nh + j] = 0.0;
1523 static real energyNoseHoover(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1527 int nh = state->nhchainlength;
1529 for (int i = 0; i < ir->opts.ngtc; i++)
1531 const double* ixi = &state->nosehoover_xi[i * nh];
1532 const double* ivxi = &state->nosehoover_vxi[i * nh];
1533 const double* iQinv = &(MassQ->Qinv[i * nh]);
1535 int nd = static_cast<int>(ir->opts.nrdf[i]);
1536 real reft = std::max<real>(ir->opts.ref_t[i], 0);
1537 real kT = BOLTZ * reft;
1541 if (inputrecNvtTrotter(ir))
1543 /* contribution from the thermal momenta of the NH chain */
1544 for (int j = 0; j < nh; j++)
1548 energy += 0.5 * gmx::square(ivxi[j]) / iQinv[j];
1549 /* contribution from the thermal variable of the NH chain */
1559 energy += ndj * ixi[j] * kT;
1563 else /* Other non Trotter temperature NH control -- no chains yet. */
1565 energy += 0.5 * BOLTZ * nd * gmx::square(ivxi[0]) / iQinv[0];
1566 energy += nd * ixi[0] * kT;
1574 /* Returns the energy from the barostat thermostat chain */
1575 static real energyPressureMTTK(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1579 int nh = state->nhchainlength;
1581 for (int i = 0; i < state->nnhpres; i++)
1583 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1584 real reft = std::max<real>(ir->opts.ref_t[0], 0.0); /* using 'System' temperature */
1585 real kT = BOLTZ * reft;
1587 for (int j = 0; j < nh; j++)
1589 double iQinv = MassQ->QPinv[i * nh + j];
1592 energy += 0.5 * gmx::square(state->nhpres_vxi[i * nh + j] / iQinv);
1593 /* contribution from the thermal variable of the NH chain */
1594 energy += state->nhpres_xi[i * nh + j] * kT;
1598 fprintf(debug, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i, j,
1599 state->nhpres_vxi[i * nh + j], state->nhpres_xi[i * nh + j]);
1607 /* Returns the energy accumulated by the V-rescale or Berendsen thermostat */
1608 static real energyVrescale(const t_inputrec* ir, const t_state* state)
1611 for (int i = 0; i < ir->opts.ngtc; i++)
1613 energy += state->therm_integral[i];
1619 real NPT_energy(const t_inputrec* ir, const t_state* state, const t_extmass* MassQ)
1623 if (ir->epc != epcNO)
1625 /* Compute the contribution of the pressure to the conserved quantity*/
1627 real vol = det(state->box);
1631 case epcPARRINELLORAHMAN:
1633 /* contribution from the pressure momenta */
1635 calcParrinelloRahmanInvMass(ir, state->box, invMass);
1636 for (int d = 0; d < DIM; d++)
1638 for (int n = 0; n <= d; n++)
1640 if (invMass[d][n] > 0)
1642 energyNPT += 0.5 * gmx::square(state->boxv[d][n]) / (invMass[d][n] * PRESFAC);
1647 /* Contribution from the PV term.
1648 * Not that with non-zero off-diagonal reference pressures,
1649 * i.e. applied shear stresses, there are additional terms.
1650 * We don't support this here, since that requires keeping
1651 * track of unwrapped box diagonal elements. This case is
1652 * excluded in integratorHasConservedEnergyQuantity().
1654 energyNPT += vol * trace(ir->ref_p) / (DIM * PRESFAC);
1658 /* contribution from the pressure momenta */
1659 energyNPT += 0.5 * gmx::square(state->veta) / MassQ->Winv;
1661 /* contribution from the PV term */
1662 energyNPT += vol * trace(ir->ref_p) / (DIM * PRESFAC);
1664 if (ir->epc == epcMTTK)
1666 /* contribution from the MTTK chain */
1667 energyNPT += energyPressureMTTK(ir, state, MassQ);
1670 case epcBERENDSEN: energyNPT += state->baros_integral; break;
1674 "Conserved energy quantity for pressure coupling is not handled. A case "
1675 "should be added with either the conserved quantity added or nothing added "
1676 "and an exclusion added to integratorHasConservedEnergyQuantity().");
1684 case etcBERENDSEN: energyNPT += energyVrescale(ir, state); break;
1685 case etcNOSEHOOVER: energyNPT += energyNoseHoover(ir, state, MassQ); break;
1687 case etcANDERSENMASSIVE:
1688 // Not supported, excluded in integratorHasConservedEnergyQuantity()
1693 "Conserved energy quantity for temperature coupling is not handled. A case "
1694 "should be added with either the conserved quantity added or nothing added and "
1695 "an exclusion added to integratorHasConservedEnergyQuantity().");
1702 static real vrescale_sumnoises(real nn, gmx::ThreeFry2x64<>* rng, gmx::NormalDistribution<real>* normalDist)
1705 * Returns the sum of nn independent gaussian noises squared
1706 * (i.e. equivalent to summing the square of the return values
1707 * of nn calls to a normal distribution).
1709 const real ndeg_tol = 0.0001;
1711 gmx::GammaDistribution<real> gammaDist(0.5 * nn, 1.0);
1713 if (nn < 2 + ndeg_tol)
1718 nn_int = gmx::roundToInt(nn);
1720 if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1723 "The v-rescale thermostat was called with a group with #DOF=%f, but for "
1724 "#DOF<3 only integer #DOF are supported",
1729 for (i = 0; i < nn_int; i++)
1731 gauss = (*normalDist)(*rng);
1737 /* Use a gamma distribution for any real nn > 2 */
1738 r = 2.0 * gammaDist(*rng);
1744 real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut, int64_t step, int64_t seed)
1747 * Generates a new value for the kinetic energy,
1748 * according to Bussi et al JCP (2007), Eq. (A7)
1749 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1750 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1751 * ndeg: number of degrees of freedom of the atoms to be thermalized
1752 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1754 real factor, rr, ekin_new;
1755 gmx::ThreeFry2x64<64> rng(seed, gmx::RandomDomain::Thermostat);
1756 gmx::NormalDistribution<real> normalDist;
1760 factor = exp(-1.0 / taut);
1767 rng.restart(step, 0);
1769 rr = normalDist(rng);
1773 * (sigma * (vrescale_sumnoises(ndeg - 1, &rng, &normalDist) + rr * rr) / ndeg - kk)
1774 + 2.0 * rr * std::sqrt(kk * sigma / ndeg * (1.0 - factor) * factor);
1779 void vrescale_tcoupl(const t_inputrec* ir, int64_t step, gmx_ekindata_t* ekind, real dt, double therm_integral[])
1781 const t_grpopts* opts;
1783 real Ek, Ek_ref1, Ek_ref, Ek_new;
1787 for (i = 0; (i < opts->ngtc); i++)
1791 Ek = trace(ekind->tcstat[i].ekinf);
1795 Ek = trace(ekind->tcstat[i].ekinh);
1798 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1800 Ek_ref1 = 0.5 * opts->ref_t[i] * BOLTZ;
1801 Ek_ref = Ek_ref1 * opts->nrdf[i];
1803 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i], opts->tau_t[i] / dt, step, ir->ld_seed);
1805 /* Analytically Ek_new>=0, but we check for rounding errors */
1808 ekind->tcstat[i].lambda = 0.0;
1812 ekind->tcstat[i].lambda = std::sqrt(Ek_new / Ek);
1815 therm_integral[i] -= Ek_new - Ek;
1819 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n", i, Ek_ref,
1820 Ek, Ek_new, ekind->tcstat[i].lambda);
1825 ekind->tcstat[i].lambda = 1.0;
1830 void rescale_velocities(const gmx_ekindata_t* ekind, const t_mdatoms* mdatoms, int start, int end, rvec v[])
1832 unsigned short *cACC, *cTC;
1839 gmx::ArrayRef<const t_grp_tcstat> tcstat = ekind->tcstat;
1843 gmx::ArrayRef<const t_grp_acc> gstat = ekind->grpstat;
1844 cACC = mdatoms->cACC;
1848 for (n = start; n < end; n++)
1858 /* Only scale the velocity component relative to the COM velocity */
1859 rvec_sub(v[n], gstat[ga].u, vrel);
1860 lg = tcstat[gt].lambda;
1861 for (d = 0; d < DIM; d++)
1863 v[n][d] = gstat[ga].u[d] + lg * vrel[d];
1870 for (n = start; n < end; n++)
1876 lg = tcstat[gt].lambda;
1877 for (d = 0; d < DIM; d++)
1885 //! Check whether we're doing simulated annealing
1886 bool doSimulatedAnnealing(const t_inputrec* ir)
1888 for (int i = 0; i < ir->opts.ngtc; i++)
1890 /* set bSimAnn if any group is being annealed */
1891 if (ir->opts.annealing[i] != eannNO)
1899 // TODO If we keep simulated annealing, make a proper module that
1900 // does not rely on changing inputrec.
1901 bool initSimulatedAnnealing(t_inputrec* ir, gmx::Update* upd)
1903 bool doSimAnnealing = doSimulatedAnnealing(ir);
1906 update_annealing_target_temp(ir, ir->init_t, upd);
1908 return doSimAnnealing;
1911 /* set target temperatures if we are annealing */
1912 void update_annealing_target_temp(t_inputrec* ir, real t, gmx::Update* upd)
1914 int i, j, n, npoints;
1915 real pert, thist = 0, x;
1917 for (i = 0; i < ir->opts.ngtc; i++)
1919 npoints = ir->opts.anneal_npoints[i];
1920 switch (ir->opts.annealing[i])
1922 case eannNO: continue;
1924 /* calculate time modulo the period */
1925 pert = ir->opts.anneal_time[i][npoints - 1];
1926 n = static_cast<int>(t / pert);
1927 thist = t - n * pert; /* modulo time */
1928 /* Make sure rounding didn't get us outside the interval */
1929 if (std::fabs(thist - pert) < GMX_REAL_EPS * 100)
1934 case eannSINGLE: thist = t; break;
1936 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",
1937 i, ir->opts.ngtc, npoints);
1939 /* We are doing annealing for this group if we got here,
1940 * and we have the (relative) time as thist.
1941 * calculate target temp */
1943 while ((j < npoints - 1) && (thist > (ir->opts.anneal_time[i][j + 1])))
1947 if (j < npoints - 1)
1949 /* Found our position between points j and j+1.
1950 * Interpolate: x is the amount from j+1, (1-x) from point j
1951 * First treat possible jumps in temperature as a special case.
1953 if ((ir->opts.anneal_time[i][j + 1] - ir->opts.anneal_time[i][j]) < GMX_REAL_EPS * 100)
1955 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][j + 1];
1959 x = ((thist - ir->opts.anneal_time[i][j])
1960 / (ir->opts.anneal_time[i][j + 1] - ir->opts.anneal_time[i][j]));
1962 x * ir->opts.anneal_temp[i][j + 1] + (1 - x) * ir->opts.anneal_temp[i][j];
1967 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][npoints - 1];
1971 upd->update_temperature_constants(*ir);
1974 void pleaseCiteCouplingAlgorithms(FILE* fplog, const t_inputrec& ir)
1976 if (EI_DYNAMICS(ir.eI))
1978 if (ir.etc == etcBERENDSEN)
1980 please_cite(fplog, "Berendsen84a");
1982 if (ir.etc == etcVRESCALE)
1984 please_cite(fplog, "Bussi2007a");
1986 // TODO this is actually an integrator, not a coupling algorithm
1989 please_cite(fplog, "Goga2012");