2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
45 #include "gromacs/domdec/domdec_struct.h"
46 #include "gromacs/gmxlib/nrnb.h"
47 #include "gromacs/math/functions.h"
48 #include "gromacs/math/invertmatrix.h"
49 #include "gromacs/math/units.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/math/vecdump.h"
52 #include "gromacs/mdlib/expanded.h"
53 #include "gromacs/mdlib/gmx_omp_nthreads.h"
54 #include "gromacs/mdlib/sim_util.h"
55 #include "gromacs/mdlib/update.h"
56 #include "gromacs/mdtypes/commrec.h"
57 #include "gromacs/mdtypes/group.h"
58 #include "gromacs/mdtypes/inputrec.h"
59 #include "gromacs/mdtypes/md_enums.h"
60 #include "gromacs/mdtypes/state.h"
61 #include "gromacs/pbcutil/boxutilities.h"
62 #include "gromacs/pbcutil/pbc.h"
63 #include "gromacs/random/gammadistribution.h"
64 #include "gromacs/random/normaldistribution.h"
65 #include "gromacs/random/tabulatednormaldistribution.h"
66 #include "gromacs/random/threefry.h"
67 #include "gromacs/random/uniformrealdistribution.h"
68 #include "gromacs/utility/cstringutil.h"
69 #include "gromacs/utility/fatalerror.h"
70 #include "gromacs/utility/smalloc.h"
72 #define NTROTTERPARTS 3
74 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
76 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
77 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
79 #define MAX_SUZUKI_YOSHIDA_NUM 5
80 #define SUZUKI_YOSHIDA_NUM 5
82 static const double sy_const_1[] = { 1. };
83 static const double sy_const_3[] = { 0.828981543588751, -0.657963087177502, 0.828981543588751 };
84 static const double sy_const_5[] = { 0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065, 0.2967324292201065 };
86 static const double* sy_const[] = {
96 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
100 {0.828981543588751,-0.657963087177502,0.828981543588751},
102 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
105 /* these integration routines are only referenced inside this file */
106 static void NHC_trotter(t_grpopts *opts, int nvar, gmx_ekindata_t *ekind, real dtfull,
107 double xi[], double vxi[], double scalefac[], real *veta, t_extmass *MassQ, gmx_bool bEkinAveVel)
110 /* general routine for both barostat and thermostat nose hoover chains */
113 double Ekin, Efac, reft, kT, nd;
115 t_grp_tcstat *tcstat;
120 int mstepsi, mstepsj;
121 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
122 int nh = opts->nhchainlength;
125 mstepsi = mstepsj = ns;
127 /* if scalefac is NULL, we are doing the NHC of the barostat */
130 if (scalefac == nullptr)
135 for (i = 0; i < nvar; i++)
138 /* make it easier to iterate by selecting
139 out the sub-array that corresponds to this T group */
145 iQinv = &(MassQ->QPinv[i*nh]);
146 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
147 reft = std::max<real>(0, opts->ref_t[0]);
148 Ekin = gmx::square(*veta)/MassQ->Winv;
152 iQinv = &(MassQ->Qinv[i*nh]);
153 tcstat = &ekind->tcstat[i];
155 reft = std::max<real>(0, opts->ref_t[i]);
158 Ekin = 2*trace(tcstat->ekinf)*tcstat->ekinscalef_nhc;
162 Ekin = 2*trace(tcstat->ekinh)*tcstat->ekinscaleh_nhc;
167 for (mi = 0; mi < mstepsi; mi++)
169 for (mj = 0; mj < mstepsj; mj++)
171 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
172 dt = sy_const[ns][mj] * dtfull / mstepsi;
174 /* compute the thermal forces */
175 GQ[0] = iQinv[0]*(Ekin - nd*kT);
177 for (j = 0; j < nh-1; j++)
181 /* we actually don't need to update here if we save the
182 state of the GQ, but it's easier to just recompute*/
183 GQ[j+1] = iQinv[j+1]*((gmx::square(ivxi[j])/iQinv[j])-kT);
191 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
192 for (j = nh-1; j > 0; j--)
194 Efac = exp(-0.125*dt*ivxi[j]);
195 ivxi[j-1] = Efac*(ivxi[j-1]*Efac + 0.25*dt*GQ[j-1]);
198 Efac = exp(-0.5*dt*ivxi[0]);
209 /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
210 able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
211 think about this a bit more . . . */
213 GQ[0] = iQinv[0]*(Ekin - nd*kT);
215 /* update thermostat positions */
216 for (j = 0; j < nh; j++)
218 ixi[j] += 0.5*dt*ivxi[j];
221 for (j = 0; j < nh-1; j++)
223 Efac = exp(-0.125*dt*ivxi[j+1]);
224 ivxi[j] = Efac*(ivxi[j]*Efac + 0.25*dt*GQ[j]);
227 GQ[j+1] = iQinv[j+1]*((gmx::square(ivxi[j])/iQinv[j])-kT);
234 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
241 static void boxv_trotter(t_inputrec *ir, real *veta, real dt, tensor box,
242 gmx_ekindata_t *ekind, tensor vir, real pcorr, t_extmass *MassQ)
249 tensor ekinmod, localpres;
251 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
252 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
255 if (ir->epct == epctSEMIISOTROPIC)
264 /* eta is in pure units. veta is in units of ps^-1. GW is in
265 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
266 taken to use only RATIOS of eta in updating the volume. */
268 /* we take the partial pressure tensors, modify the
269 kinetic energy tensor, and recovert to pressure */
271 if (ir->opts.nrdf[0] == 0)
273 gmx_fatal(FARGS, "Barostat is coupled to a T-group with no degrees of freedom\n");
275 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
276 alpha = 1.0 + DIM/(static_cast<double>(ir->opts.nrdf[0]));
277 alpha *= ekind->tcstat[0].ekinscalef_nhc;
278 msmul(ekind->ekin, alpha, ekinmod);
279 /* for now, we use Elr = 0, because if you want to get it right, you
280 really should be using PME. Maybe print a warning? */
282 pscal = calc_pres(ir->ePBC, nwall, box, ekinmod, vir, localpres)+pcorr;
285 GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
291 * This file implements temperature and pressure coupling algorithms:
292 * For now only the Weak coupling and the modified weak coupling.
294 * Furthermore computation of pressure and temperature is done here
298 real calc_pres(int ePBC, int nwall, matrix box, tensor ekin, tensor vir,
304 if (ePBC == epbcNONE || (ePBC == epbcXY && nwall != 2))
310 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
311 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
315 fac = PRESFAC*2.0/det(box);
316 for (n = 0; (n < DIM); n++)
318 for (m = 0; (m < DIM); m++)
320 pres[n][m] = (ekin[n][m] - vir[n][m])*fac;
326 pr_rvecs(debug, 0, "PC: pres", pres, DIM);
327 pr_rvecs(debug, 0, "PC: ekin", ekin, DIM);
328 pr_rvecs(debug, 0, "PC: vir ", vir, DIM);
329 pr_rvecs(debug, 0, "PC: box ", box, DIM);
332 return trace(pres)/DIM;
335 real calc_temp(real ekin, real nrdf)
339 return (2.0*ekin)/(nrdf*BOLTZ);
347 /*! \brief Sets 1/mass for Parrinello-Rahman in wInv; NOTE: PRESFAC is not included, so not in GROMACS units! */
348 static void calcParrinelloRahmanInvMass(const t_inputrec *ir, const matrix box,
353 /* TODO: See if we can make the mass independent of the box size */
354 maxBoxLength = std::max(box[XX][XX], box[YY][YY]);
355 maxBoxLength = std::max(maxBoxLength, box[ZZ][ZZ]);
357 for (int d = 0; d < DIM; d++)
359 for (int n = 0; n < DIM; n++)
361 wInv[d][n] = (4*M_PI*M_PI*ir->compress[d][n])/(3*ir->tau_p*ir->tau_p*maxBoxLength);
366 void parrinellorahman_pcoupl(FILE *fplog, int64_t step,
367 const t_inputrec *ir, real dt, const tensor pres,
368 tensor box, tensor box_rel, tensor boxv,
369 tensor M, matrix mu, gmx_bool bFirstStep)
371 /* This doesn't do any coordinate updating. It just
372 * integrates the box vector equations from the calculated
373 * acceleration due to pressure difference. We also compute
374 * the tensor M which is used in update to couple the particle
375 * coordinates to the box vectors.
377 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
380 * M_nk = (h') * (h' * h + h' h) * h
382 * with the dots denoting time derivatives and h is the transformation from
383 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
384 * This also goes for the pressure and M tensors - they are transposed relative
385 * to ours. Our equation thus becomes:
388 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
390 * where b is the gromacs box matrix.
391 * Our box accelerations are given by
393 * b = vol/W inv(box') * (P-ref_P) (=h')
396 real vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
397 real atot, arel, change;
398 tensor invbox, pdiff, t1, t2;
400 gmx::invertBoxMatrix(box, invbox);
404 /* Note that PRESFAC does not occur here.
405 * The pressure and compressibility always occur as a product,
406 * therefore the pressure unit drops out.
409 calcParrinelloRahmanInvMass(ir, box, winv);
411 m_sub(pres, ir->ref_p, pdiff);
413 if (ir->epct == epctSURFACETENSION)
415 /* Unlike Berendsen coupling it might not be trivial to include a z
416 * pressure correction here? On the other hand we don't scale the
417 * box momentarily, but change accelerations, so it might not be crucial.
419 real xy_pressure = 0.5*(pres[XX][XX]+pres[YY][YY]);
420 for (int d = 0; d < ZZ; d++)
422 pdiff[d][d] = (xy_pressure-(pres[ZZ][ZZ]-ir->ref_p[d][d]/box[d][d]));
426 tmmul(invbox, pdiff, t1);
427 /* Move the off-diagonal elements of the 'force' to one side to ensure
428 * that we obey the box constraints.
430 for (int d = 0; d < DIM; d++)
432 for (int n = 0; n < d; n++)
434 t1[d][n] += t1[n][d];
441 case epctANISOTROPIC:
442 for (int d = 0; d < DIM; d++)
444 for (int n = 0; n <= d; n++)
446 t1[d][n] *= winv[d][n]*vol;
451 /* calculate total volume acceleration */
452 atot = box[XX][XX]*box[YY][YY]*t1[ZZ][ZZ]+
453 box[XX][XX]*t1[YY][YY]*box[ZZ][ZZ]+
454 t1[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
456 /* set all RELATIVE box accelerations equal, and maintain total V
458 for (int d = 0; d < DIM; d++)
460 for (int n = 0; n <= d; n++)
462 t1[d][n] = winv[0][0]*vol*arel*box[d][n];
466 case epctSEMIISOTROPIC:
467 case epctSURFACETENSION:
468 /* Note the correction to pdiff above for surftens. coupling */
470 /* calculate total XY volume acceleration */
471 atot = box[XX][XX]*t1[YY][YY]+t1[XX][XX]*box[YY][YY];
472 arel = atot/(2*box[XX][XX]*box[YY][YY]);
473 /* set RELATIVE XY box accelerations equal, and maintain total V
474 * change speed. Dont change the third box vector accelerations */
475 for (int d = 0; d < ZZ; d++)
477 for (int n = 0; n <= d; n++)
479 t1[d][n] = winv[d][n]*vol*arel*box[d][n];
482 for (int n = 0; n < DIM; n++)
484 t1[ZZ][n] *= winv[ZZ][n]*vol;
488 gmx_fatal(FARGS, "Parrinello-Rahman pressure coupling type %s "
489 "not supported yet\n", EPCOUPLTYPETYPE(ir->epct));
493 for (int d = 0; d < DIM; d++)
495 for (int n = 0; n <= d; n++)
497 boxv[d][n] += dt*t1[d][n];
499 /* We do NOT update the box vectors themselves here, since
500 * we need them for shifting later. It is instead done last
501 * in the update() routine.
504 /* Calculate the change relative to diagonal elements-
505 since it's perfectly ok for the off-diagonal ones to
506 be zero it doesn't make sense to check the change relative
510 change = std::fabs(dt*boxv[d][n]/box[d][d]);
512 if (change > maxchange)
519 if (maxchange > 0.01 && fplog)
523 "\nStep %s Warning: Pressure scaling more than 1%%. "
524 "This may mean your system\n is not yet equilibrated. "
525 "Use of Parrinello-Rahman pressure coupling during\n"
526 "equilibration can lead to simulation instability, "
527 "and is discouraged.\n",
528 gmx_step_str(step, buf));
532 preserve_box_shape(ir, box_rel, boxv);
534 mtmul(boxv, box, t1); /* t1=boxv * b' */
535 mmul(invbox, t1, t2);
536 mtmul(t2, invbox, M);
538 /* Determine the scaling matrix mu for the coordinates */
539 for (int d = 0; d < DIM; d++)
541 for (int n = 0; n <= d; n++)
543 t1[d][n] = box[d][n] + dt*boxv[d][n];
546 preserve_box_shape(ir, box_rel, t1);
547 /* t1 is the box at t+dt, determine mu as the relative change */
548 mmul_ur0(invbox, t1, mu);
551 void berendsen_pcoupl(FILE *fplog, int64_t step,
552 const t_inputrec *ir, real dt,
553 const tensor pres, const matrix box,
554 const matrix force_vir, const matrix constraint_vir,
555 matrix mu, double *baros_integral)
558 real scalar_pressure, xy_pressure, p_corr_z;
562 * Calculate the scaling matrix mu
566 for (d = 0; d < DIM; d++)
568 scalar_pressure += pres[d][d]/DIM;
571 xy_pressure += pres[d][d]/(DIM-1);
574 /* Pressure is now in bar, everywhere. */
575 #define factor(d, m) (ir->compress[d][m]*dt/ir->tau_p)
577 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
578 * necessary for triclinic scaling
584 for (d = 0; d < DIM; d++)
586 mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
589 case epctSEMIISOTROPIC:
590 for (d = 0; d < ZZ; d++)
592 mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
595 1.0 - factor(ZZ, ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
597 case epctANISOTROPIC:
598 for (d = 0; d < DIM; d++)
600 for (n = 0; n < DIM; n++)
602 mu[d][n] = (d == n ? 1.0 : 0.0)
603 -factor(d, n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
607 case epctSURFACETENSION:
608 /* ir->ref_p[0/1] is the reference surface-tension times *
609 * the number of surfaces */
610 if (ir->compress[ZZ][ZZ])
612 p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
616 /* when the compressibity is zero, set the pressure correction *
617 * in the z-direction to zero to get the correct surface tension */
620 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
621 for (d = 0; d < DIM-1; d++)
623 mu[d][d] = 1.0 + factor(d, d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
624 - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
628 gmx_fatal(FARGS, "Berendsen pressure coupling type %s not supported yet\n",
629 EPCOUPLTYPETYPE(ir->epct));
631 /* To fullfill the orientation restrictions on triclinic boxes
632 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
633 * the other elements of mu to first order.
635 mu[YY][XX] += mu[XX][YY];
636 mu[ZZ][XX] += mu[XX][ZZ];
637 mu[ZZ][YY] += mu[YY][ZZ];
642 /* Keep track of the work the barostat applies on the system.
643 * Without constraints force_vir tells us how Epot changes when scaling.
644 * With constraints constraint_vir gives us the constraint contribution
645 * to both Epot and Ekin. Although we are not scaling velocities, scaling
646 * the coordinates leads to scaling of distances involved in constraints.
647 * This in turn changes the angular momentum (even if the constrained
648 * distances are corrected at the next step). The kinetic component
649 * of the constraint virial captures the angular momentum change.
651 for (int d = 0; d < DIM; d++)
653 for (int n = 0; n <= d; n++)
655 *baros_integral -= 2*(mu[d][n] - (n == d ? 1 : 0))*(force_vir[d][n] + constraint_vir[d][n]);
661 pr_rvecs(debug, 0, "PC: pres ", pres, 3);
662 pr_rvecs(debug, 0, "PC: mu ", mu, 3);
665 if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 ||
666 mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01 ||
667 mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
670 sprintf(buf, "\nStep %s Warning: pressure scaling more than 1%%, "
672 gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
675 fprintf(fplog, "%s", buf);
677 fprintf(stderr, "%s", buf);
681 void berendsen_pscale(const t_inputrec *ir, const matrix mu,
682 matrix box, matrix box_rel,
683 int start, int nr_atoms,
684 rvec x[], const unsigned short cFREEZE[],
687 ivec *nFreeze = ir->opts.nFreeze;
689 int nthreads gmx_unused;
691 #ifndef __clang_analyzer__
692 nthreads = gmx_omp_nthreads_get(emntUpdate);
695 /* Scale the positions */
696 #pragma omp parallel for num_threads(nthreads) schedule(static)
697 for (n = start; n < start+nr_atoms; n++)
699 // Trivial OpenMP region that does not throw
702 if (cFREEZE == nullptr)
713 x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
717 x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
721 x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
724 /* compute final boxlengths */
725 for (d = 0; d < DIM; d++)
727 box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
728 box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
729 box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
732 preserve_box_shape(ir, box_rel, box);
734 /* (un)shifting should NOT be done after this,
735 * since the box vectors might have changed
737 inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
740 void berendsen_tcoupl(const t_inputrec *ir, const gmx_ekindata_t *ekind, real dt,
741 std::vector<double> &therm_integral)
743 const t_grpopts *opts = &ir->opts;
745 for (int i = 0; (i < opts->ngtc); i++)
751 Ek = trace(ekind->tcstat[i].ekinf);
752 T = ekind->tcstat[i].T;
756 Ek = trace(ekind->tcstat[i].ekinh);
757 T = ekind->tcstat[i].Th;
760 if ((opts->tau_t[i] > 0) && (T > 0.0))
762 real reft = std::max<real>(0, opts->ref_t[i]);
763 real lll = std::sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
764 ekind->tcstat[i].lambda = std::max<real>(std::min<real>(lll, 1.25), 0.8);
768 ekind->tcstat[i].lambda = 1.0;
771 /* Keep track of the amount of energy we are adding to the system */
772 therm_integral[i] -= (gmx::square(ekind->tcstat[i].lambda) - 1)*Ek;
776 fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n",
777 i, T, ekind->tcstat[i].lambda);
782 void andersen_tcoupl(t_inputrec *ir, int64_t step,
783 const t_commrec *cr, const t_mdatoms *md, t_state *state, real rate, const gmx_bool *randomize, const real *boltzfac)
785 const int *gatindex = (DOMAINDECOMP(cr) ? cr->dd->globalAtomIndices.data() : nullptr);
788 gmx::ThreeFry2x64<0> rng(ir->andersen_seed, gmx::RandomDomain::Thermostat);
789 gmx::UniformRealDistribution<real> uniformDist;
790 gmx::TabulatedNormalDistribution<real, 14> normalDist;
792 /* randomize the velocities of the selected particles */
794 for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */
796 int ng = gatindex ? gatindex[i] : i;
799 rng.restart(step, ng);
803 gc = md->cTC[i]; /* assign the atom to a temperature group if there are more than one */
807 if (ir->etc == etcANDERSENMASSIVE)
809 /* Randomize particle always */
814 /* Randomize particle probabilistically */
816 bRandomize = uniformDist(rng) < rate;
823 scal = std::sqrt(boltzfac[gc]*md->invmass[i]);
827 for (d = 0; d < DIM; d++)
829 state->v[i][d] = scal*normalDist(rng);
837 void nosehoover_tcoupl(t_grpopts *opts, gmx_ekindata_t *ekind, real dt,
838 double xi[], double vxi[], t_extmass *MassQ)
843 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
845 for (i = 0; (i < opts->ngtc); i++)
847 reft = std::max<real>(0, opts->ref_t[i]);
849 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
850 xi[i] += dt*(oldvxi + vxi[i])*0.5;
854 void trotter_update(t_inputrec *ir, int64_t step, gmx_ekindata_t *ekind,
855 gmx_enerdata_t *enerd, t_state *state,
856 tensor vir, t_mdatoms *md,
857 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
860 int n, i, d, ngtc, gc = 0, t;
861 t_grp_tcstat *tcstat;
865 double *scalefac, dtc;
867 rvec sumv = {0, 0, 0};
870 if (trotter_seqno <= ettTSEQ2)
872 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
873 is actually the last half step from the previous step. Thus the first half step
874 actually corresponds to the n-1 step*/
882 bCouple = (ir->nsttcouple == 1 ||
883 do_per_step(step_eff+ir->nsttcouple, ir->nsttcouple));
885 trotter_seq = trotter_seqlist[trotter_seqno];
887 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
891 dtc = ir->nsttcouple*ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
892 opts = &(ir->opts); /* just for ease of referencing */
895 snew(scalefac, opts->ngtc);
896 for (i = 0; i < ngtc; i++)
900 /* execute the series of trotter updates specified in the trotterpart array */
902 for (i = 0; i < NTROTTERPARTS; i++)
904 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
905 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
914 switch (trotter_seq[i])
918 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir,
919 enerd->term[F_PDISPCORR], MassQ);
923 NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi.data(),
924 state->nhpres_vxi.data(), nullptr, &(state->veta), MassQ, FALSE);
928 NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi.data(),
929 state->nosehoover_vxi.data(), scalefac, nullptr, MassQ, (ir->eI == eiVV));
930 /* need to rescale the kinetic energies and velocities here. Could
931 scale the velocities later, but we need them scaled in order to
932 produce the correct outputs, so we'll scale them here. */
934 for (t = 0; t < ngtc; t++)
936 tcstat = &ekind->tcstat[t];
937 tcstat->vscale_nhc = scalefac[t];
938 tcstat->ekinscaleh_nhc *= (scalefac[t]*scalefac[t]);
939 tcstat->ekinscalef_nhc *= (scalefac[t]*scalefac[t]);
941 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
942 /* but do we actually need the total? */
944 /* modify the velocities as well */
945 for (n = 0; n < md->homenr; n++)
947 if (md->cTC) /* does this conditional need to be here? is this always true?*/
951 for (d = 0; d < DIM; d++)
953 state->v[n][d] *= scalefac[gc];
958 for (d = 0; d < DIM; d++)
960 sumv[d] += (state->v[n][d])/md->invmass[n];
969 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
974 extern void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit)
976 int n, i, j, d, ngtc, nh;
978 real reft, kT, ndj, nd;
980 opts = &(ir->opts); /* just for ease of referencing */
981 ngtc = ir->opts.ngtc;
982 nh = state->nhchainlength;
988 snew(MassQ->Qinv, ngtc);
990 for (i = 0; (i < ngtc); i++)
992 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
994 MassQ->Qinv[i] = 1.0/(gmx::square(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
998 MassQ->Qinv[i] = 0.0;
1002 else if (EI_VV(ir->eI))
1004 /* Set pressure variables */
1008 if (state->vol0 == 0)
1010 state->vol0 = det(state->box);
1011 /* because we start by defining a fixed
1012 compressibility, we need the volume at this
1013 compressibility to solve the problem. */
1017 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1018 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1019 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*gmx::square(ir->tau_p/M_2PI));
1020 /* An alternate mass definition, from Tuckerman et al. */
1021 /* MassQ->Winv = 1.0/(gmx::square(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1022 for (d = 0; d < DIM; d++)
1024 for (n = 0; n < DIM; n++)
1026 MassQ->Winvm[d][n] = PRESFAC*ir->compress[d][n]/(state->vol0*gmx::square(ir->tau_p/M_2PI));
1027 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1028 before using MTTK for anisotropic states.*/
1031 /* Allocate space for thermostat variables */
1034 snew(MassQ->Qinv, ngtc*nh);
1037 /* now, set temperature variables */
1038 for (i = 0; i < ngtc; i++)
1040 if (opts->tau_t[i] > 0 && opts->ref_t[i] > 0 && opts->nrdf[i] > 0)
1042 reft = std::max<real>(0, opts->ref_t[i]);
1045 for (j = 0; j < nh; j++)
1055 MassQ->Qinv[i*nh+j] = 1.0/(gmx::square(opts->tau_t[i]/M_2PI)*ndj*kT);
1060 for (j = 0; j < nh; j++)
1062 MassQ->Qinv[i*nh+j] = 0.0;
1069 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
1071 int i, j, nnhpres, nh;
1073 real bmass, qmass, reft, kT;
1076 opts = &(ir->opts); /* just for ease of referencing */
1077 nnhpres = state->nnhpres;
1078 nh = state->nhchainlength;
1080 if (EI_VV(ir->eI) && (ir->epc == epcMTTK) && (ir->etc != etcNOSEHOOVER))
1082 gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
1085 init_npt_masses(ir, state, MassQ, TRUE);
1087 /* first, initialize clear all the trotter calls */
1088 snew(trotter_seq, ettTSEQMAX);
1089 for (i = 0; i < ettTSEQMAX; i++)
1091 snew(trotter_seq[i], NTROTTERPARTS);
1092 for (j = 0; j < NTROTTERPARTS; j++)
1094 trotter_seq[i][j] = etrtNONE;
1096 trotter_seq[i][0] = etrtSKIPALL;
1101 /* no trotter calls, so we never use the values in the array.
1102 * We access them (so we need to define them, but ignore
1108 /* compute the kinetic energy by using the half step velocities or
1109 * the kinetic energies, depending on the order of the trotter calls */
1113 if (inputrecNptTrotter(ir))
1115 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1116 We start with the initial one. */
1117 /* first, a round that estimates veta. */
1118 trotter_seq[0][0] = etrtBAROV;
1120 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1122 /* The first half trotter update */
1123 trotter_seq[2][0] = etrtBAROV;
1124 trotter_seq[2][1] = etrtNHC;
1125 trotter_seq[2][2] = etrtBARONHC;
1127 /* The second half trotter update */
1128 trotter_seq[3][0] = etrtBARONHC;
1129 trotter_seq[3][1] = etrtNHC;
1130 trotter_seq[3][2] = etrtBAROV;
1132 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1135 else if (inputrecNvtTrotter(ir))
1137 /* This is the easy version - there are only two calls, both the same.
1138 Otherwise, even easier -- no calls */
1139 trotter_seq[2][0] = etrtNHC;
1140 trotter_seq[3][0] = etrtNHC;
1142 else if (inputrecNphTrotter(ir))
1144 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1145 We start with the initial one. */
1146 /* first, a round that estimates veta. */
1147 trotter_seq[0][0] = etrtBAROV;
1149 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1151 /* The first half trotter update */
1152 trotter_seq[2][0] = etrtBAROV;
1153 trotter_seq[2][1] = etrtBARONHC;
1155 /* The second half trotter update */
1156 trotter_seq[3][0] = etrtBARONHC;
1157 trotter_seq[3][1] = etrtBAROV;
1159 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1162 else if (ir->eI == eiVVAK)
1164 if (inputrecNptTrotter(ir))
1166 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1167 We start with the initial one. */
1168 /* first, a round that estimates veta. */
1169 trotter_seq[0][0] = etrtBAROV;
1171 /* The first half trotter update, part 1 -- double update, because it commutes */
1172 trotter_seq[1][0] = etrtNHC;
1174 /* The first half trotter update, part 2 */
1175 trotter_seq[2][0] = etrtBAROV;
1176 trotter_seq[2][1] = etrtBARONHC;
1178 /* The second half trotter update, part 1 */
1179 trotter_seq[3][0] = etrtBARONHC;
1180 trotter_seq[3][1] = etrtBAROV;
1182 /* The second half trotter update */
1183 trotter_seq[4][0] = etrtNHC;
1185 else if (inputrecNvtTrotter(ir))
1187 /* This is the easy version - there is only one call, both the same.
1188 Otherwise, even easier -- no calls */
1189 trotter_seq[1][0] = etrtNHC;
1190 trotter_seq[4][0] = etrtNHC;
1192 else if (inputrecNphTrotter(ir))
1194 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1195 We start with the initial one. */
1196 /* first, a round that estimates veta. */
1197 trotter_seq[0][0] = etrtBAROV;
1199 /* The first half trotter update, part 1 -- leave zero */
1200 trotter_seq[1][0] = etrtNHC;
1202 /* The first half trotter update, part 2 */
1203 trotter_seq[2][0] = etrtBAROV;
1204 trotter_seq[2][1] = etrtBARONHC;
1206 /* The second half trotter update, part 1 */
1207 trotter_seq[3][0] = etrtBARONHC;
1208 trotter_seq[3][1] = etrtBAROV;
1210 /* The second half trotter update -- blank for now */
1218 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
1221 snew(MassQ->QPinv, nnhpres*opts->nhchainlength);
1223 /* barostat temperature */
1224 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1226 reft = std::max<real>(0, opts->ref_t[0]);
1228 for (i = 0; i < nnhpres; i++)
1230 for (j = 0; j < nh; j++)
1240 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(gmx::square(opts->tau_t[0]/M_2PI)*qmass*kT);
1246 for (i = 0; i < nnhpres; i++)
1248 for (j = 0; j < nh; j++)
1250 MassQ->QPinv[i*nh+j] = 0.0;
1257 static real energyNoseHoover(const t_inputrec *ir, const t_state *state, const t_extmass *MassQ)
1261 int nh = state->nhchainlength;
1263 for (int i = 0; i < ir->opts.ngtc; i++)
1265 const double *ixi = &state->nosehoover_xi[i*nh];
1266 const double *ivxi = &state->nosehoover_vxi[i*nh];
1267 const double *iQinv = &(MassQ->Qinv[i*nh]);
1269 int nd = ir->opts.nrdf[i];
1270 real reft = std::max<real>(ir->opts.ref_t[i], 0);
1271 real kT = BOLTZ * reft;
1275 if (inputrecNvtTrotter(ir))
1277 /* contribution from the thermal momenta of the NH chain */
1278 for (int j = 0; j < nh; j++)
1282 energy += 0.5*gmx::square(ivxi[j])/iQinv[j];
1283 /* contribution from the thermal variable of the NH chain */
1293 energy += ndj*ixi[j]*kT;
1297 else /* Other non Trotter temperature NH control -- no chains yet. */
1299 energy += 0.5*BOLTZ*nd*gmx::square(ivxi[0])/iQinv[0];
1300 energy += nd*ixi[0]*kT;
1308 /* Returns the energy from the barostat thermostat chain */
1309 static real energyPressureMTTK(const t_inputrec *ir, const t_state *state, const t_extmass *MassQ)
1313 int nh = state->nhchainlength;
1315 for (int i = 0; i < state->nnhpres; i++)
1317 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1318 real reft = std::max<real>(ir->opts.ref_t[0], 0.0); /* using 'System' temperature */
1319 real kT = BOLTZ * reft;
1321 for (int j = 0; j < nh; j++)
1323 double iQinv = MassQ->QPinv[i*nh + j];
1326 energy += 0.5*gmx::square(state->nhpres_vxi[i*nh + j]/iQinv);
1327 /* contribution from the thermal variable of the NH chain */
1328 energy += state->nhpres_xi[i*nh + j]*kT;
1332 fprintf(debug, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i, j, state->nhpres_vxi[i*nh + j], state->nhpres_xi[i*nh + j]);
1340 /* Returns the energy accumulated by the V-rescale or Berendsen thermostat */
1341 static real energyVrescale(const t_inputrec *ir, const t_state *state)
1344 for (int i = 0; i < ir->opts.ngtc; i++)
1346 energy += state->therm_integral[i];
1352 real NPT_energy(const t_inputrec *ir, const t_state *state, const t_extmass *MassQ)
1356 if (ir->epc != epcNO)
1358 /* Compute the contribution of the pressure to the conserved quantity*/
1360 real vol = det(state->box);
1364 case epcPARRINELLORAHMAN:
1366 /* contribution from the pressure momenta */
1368 calcParrinelloRahmanInvMass(ir, state->box, invMass);
1369 for (int d = 0; d < DIM; d++)
1371 for (int n = 0; n <= d; n++)
1373 if (invMass[d][n] > 0)
1375 energyNPT += 0.5*gmx::square(state->boxv[d][n])/(invMass[d][n]*PRESFAC);
1380 /* Contribution from the PV term.
1381 * Not that with non-zero off-diagonal reference pressures,
1382 * i.e. applied shear stresses, there are additional terms.
1383 * We don't support this here, since that requires keeping
1384 * track of unwrapped box diagonal elements. This case is
1385 * excluded in integratorHasConservedEnergyQuantity().
1387 energyNPT += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1391 /* contribution from the pressure momenta */
1392 energyNPT += 0.5*gmx::square(state->veta)/MassQ->Winv;
1394 /* contribution from the PV term */
1395 energyNPT += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1397 if (ir->epc == epcMTTK)
1399 /* contribution from the MTTK chain */
1400 energyNPT += energyPressureMTTK(ir, state, MassQ);
1404 energyNPT += state->baros_integral;
1407 GMX_RELEASE_ASSERT(false, "Conserved energy quantity for pressure coupling is not handled. A case should be added with either the conserved quantity added or nothing added and an exclusion added to integratorHasConservedEnergyQuantity().");
1417 energyNPT += energyVrescale(ir, state);
1420 energyNPT += energyNoseHoover(ir, state, MassQ);
1423 case etcANDERSENMASSIVE:
1424 // Not supported, excluded in integratorHasConservedEnergyQuantity()
1427 GMX_RELEASE_ASSERT(false, "Conserved energy quantity for temperature coupling is not handled. A case should be added with either the conserved quantity added or nothing added and an exclusion added to integratorHasConservedEnergyQuantity().");
1434 static real vrescale_sumnoises(real nn,
1435 gmx::ThreeFry2x64<> *rng,
1436 gmx::NormalDistribution<real> *normalDist)
1439 * Returns the sum of nn independent gaussian noises squared
1440 * (i.e. equivalent to summing the square of the return values
1441 * of nn calls to a normal distribution).
1443 const real ndeg_tol = 0.0001;
1445 gmx::GammaDistribution<real> gammaDist(0.5*nn, 1.0);
1447 if (nn < 2 + ndeg_tol)
1452 nn_int = static_cast<int>(nn + 0.5);
1454 if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1456 gmx_fatal(FARGS, "The v-rescale thermostat was called with a group with #DOF=%f, but for #DOF<3 only integer #DOF are supported", nn + 1);
1460 for (i = 0; i < nn_int; i++)
1462 gauss = (*normalDist)(*rng);
1468 /* Use a gamma distribution for any real nn > 2 */
1469 r = 2.0*gammaDist(*rng);
1475 static real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut,
1476 int64_t step, int64_t seed)
1479 * Generates a new value for the kinetic energy,
1480 * according to Bussi et al JCP (2007), Eq. (A7)
1481 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1482 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1483 * ndeg: number of degrees of freedom of the atoms to be thermalized
1484 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1486 real factor, rr, ekin_new;
1487 gmx::ThreeFry2x64<64> rng(seed, gmx::RandomDomain::Thermostat);
1488 gmx::NormalDistribution<real> normalDist;
1492 factor = exp(-1.0/taut);
1499 rng.restart(step, 0);
1501 rr = normalDist(rng);
1505 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, &rng, &normalDist) + rr*rr)/ndeg - kk) +
1506 2.0*rr*std::sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1511 void vrescale_tcoupl(t_inputrec *ir, int64_t step,
1512 gmx_ekindata_t *ekind, real dt,
1513 double therm_integral[])
1517 real Ek, Ek_ref1, Ek_ref, Ek_new;
1521 for (i = 0; (i < opts->ngtc); i++)
1525 Ek = trace(ekind->tcstat[i].ekinf);
1529 Ek = trace(ekind->tcstat[i].ekinh);
1532 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1534 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1535 Ek_ref = Ek_ref1*opts->nrdf[i];
1537 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i],
1541 /* Analytically Ek_new>=0, but we check for rounding errors */
1544 ekind->tcstat[i].lambda = 0.0;
1548 ekind->tcstat[i].lambda = std::sqrt(Ek_new/Ek);
1551 therm_integral[i] -= Ek_new - Ek;
1555 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1556 i, Ek_ref, Ek, Ek_new, ekind->tcstat[i].lambda);
1561 ekind->tcstat[i].lambda = 1.0;
1566 void rescale_velocities(gmx_ekindata_t *ekind, t_mdatoms *mdatoms,
1567 int start, int end, rvec v[])
1570 t_grp_tcstat *tcstat;
1571 unsigned short *cACC, *cTC;
1576 tcstat = ekind->tcstat;
1581 gstat = ekind->grpstat;
1582 cACC = mdatoms->cACC;
1586 for (n = start; n < end; n++)
1596 /* Only scale the velocity component relative to the COM velocity */
1597 rvec_sub(v[n], gstat[ga].u, vrel);
1598 lg = tcstat[gt].lambda;
1599 for (d = 0; d < DIM; d++)
1601 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1608 for (n = start; n < end; n++)
1614 lg = tcstat[gt].lambda;
1615 for (d = 0; d < DIM; d++)
1624 /* set target temperatures if we are annealing */
1625 void update_annealing_target_temp(t_inputrec *ir, real t, gmx_update_t *upd)
1627 int i, j, n, npoints;
1628 real pert, thist = 0, x;
1630 for (i = 0; i < ir->opts.ngtc; i++)
1632 npoints = ir->opts.anneal_npoints[i];
1633 switch (ir->opts.annealing[i])
1638 /* calculate time modulo the period */
1639 pert = ir->opts.anneal_time[i][npoints-1];
1640 n = static_cast<int>(t / pert);
1641 thist = t - n*pert; /* modulo time */
1642 /* Make sure rounding didn't get us outside the interval */
1643 if (std::fabs(thist-pert) < GMX_REAL_EPS*100)
1652 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)", i, ir->opts.ngtc, npoints);
1654 /* We are doing annealing for this group if we got here,
1655 * and we have the (relative) time as thist.
1656 * calculate target temp */
1658 while ((j < npoints-1) && (thist > (ir->opts.anneal_time[i][j+1])))
1664 /* Found our position between points j and j+1.
1665 * Interpolate: x is the amount from j+1, (1-x) from point j
1666 * First treat possible jumps in temperature as a special case.
1668 if ((ir->opts.anneal_time[i][j+1]-ir->opts.anneal_time[i][j]) < GMX_REAL_EPS*100)
1670 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][j+1];
1674 x = ((thist-ir->opts.anneal_time[i][j])/
1675 (ir->opts.anneal_time[i][j+1]-ir->opts.anneal_time[i][j]));
1676 ir->opts.ref_t[i] = x*ir->opts.anneal_temp[i][j+1]+(1-x)*ir->opts.anneal_temp[i][j];
1681 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][npoints-1];
1685 update_temperature_constants(upd, ir);