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, 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.
43 #include "gromacs/domdec/domdec_struct.h"
44 #include "gromacs/gmxlib/nrnb.h"
45 #include "gromacs/math/functions.h"
46 #include "gromacs/math/invertmatrix.h"
47 #include "gromacs/math/units.h"
48 #include "gromacs/math/vec.h"
49 #include "gromacs/math/vecdump.h"
50 #include "gromacs/mdlib/gmx_omp_nthreads.h"
51 #include "gromacs/mdlib/mdrun.h"
52 #include "gromacs/mdlib/sim_util.h"
53 #include "gromacs/mdlib/update.h"
54 #include "gromacs/mdtypes/commrec.h"
55 #include "gromacs/mdtypes/group.h"
56 #include "gromacs/mdtypes/inputrec.h"
57 #include "gromacs/mdtypes/md_enums.h"
58 #include "gromacs/pbcutil/boxutilities.h"
59 #include "gromacs/pbcutil/pbc.h"
60 #include "gromacs/random/gammadistribution.h"
61 #include "gromacs/random/normaldistribution.h"
62 #include "gromacs/random/tabulatednormaldistribution.h"
63 #include "gromacs/random/threefry.h"
64 #include "gromacs/random/uniformrealdistribution.h"
65 #include "gromacs/trajectory/energy.h"
66 #include "gromacs/utility/cstringutil.h"
67 #include "gromacs/utility/fatalerror.h"
68 #include "gromacs/utility/smalloc.h"
70 #define NTROTTERPARTS 3
72 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
74 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
75 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
77 #define MAX_SUZUKI_YOSHIDA_NUM 5
78 #define SUZUKI_YOSHIDA_NUM 5
80 static const double sy_const_1[] = { 1. };
81 static const double sy_const_3[] = { 0.828981543588751, -0.657963087177502, 0.828981543588751 };
82 static const double sy_const_5[] = { 0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065, 0.2967324292201065 };
84 static const double* sy_const[] = {
94 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
98 {0.828981543588751,-0.657963087177502,0.828981543588751},
100 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
103 /* these integration routines are only referenced inside this file */
104 static void NHC_trotter(t_grpopts *opts, int nvar, gmx_ekindata_t *ekind, real dtfull,
105 double xi[], double vxi[], double scalefac[], real *veta, t_extmass *MassQ, gmx_bool bEkinAveVel)
108 /* general routine for both barostat and thermostat nose hoover chains */
111 double Ekin, Efac, reft, kT, nd;
113 t_grp_tcstat *tcstat;
118 int mstepsi, mstepsj;
119 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
120 int nh = opts->nhchainlength;
123 mstepsi = mstepsj = ns;
125 /* if scalefac is NULL, we are doing the NHC of the barostat */
128 if (scalefac == NULL)
133 for (i = 0; i < nvar; i++)
136 /* make it easier to iterate by selecting
137 out the sub-array that corresponds to this T group */
143 iQinv = &(MassQ->QPinv[i*nh]);
144 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
145 reft = std::max<real>(0, opts->ref_t[0]);
146 Ekin = gmx::square(*veta)/MassQ->Winv;
150 iQinv = &(MassQ->Qinv[i*nh]);
151 tcstat = &ekind->tcstat[i];
153 reft = std::max<real>(0, opts->ref_t[i]);
156 Ekin = 2*trace(tcstat->ekinf)*tcstat->ekinscalef_nhc;
160 Ekin = 2*trace(tcstat->ekinh)*tcstat->ekinscaleh_nhc;
165 for (mi = 0; mi < mstepsi; mi++)
167 for (mj = 0; mj < mstepsj; mj++)
169 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
170 dt = sy_const[ns][mj] * dtfull / mstepsi;
172 /* compute the thermal forces */
173 GQ[0] = iQinv[0]*(Ekin - nd*kT);
175 for (j = 0; j < nh-1; j++)
179 /* we actually don't need to update here if we save the
180 state of the GQ, but it's easier to just recompute*/
181 GQ[j+1] = iQinv[j+1]*((gmx::square(ivxi[j])/iQinv[j])-kT);
189 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
190 for (j = nh-1; j > 0; j--)
192 Efac = exp(-0.125*dt*ivxi[j]);
193 ivxi[j-1] = Efac*(ivxi[j-1]*Efac + 0.25*dt*GQ[j-1]);
196 Efac = exp(-0.5*dt*ivxi[0]);
207 /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
208 able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
209 think about this a bit more . . . */
211 GQ[0] = iQinv[0]*(Ekin - nd*kT);
213 /* update thermostat positions */
214 for (j = 0; j < nh; j++)
216 ixi[j] += 0.5*dt*ivxi[j];
219 for (j = 0; j < nh-1; j++)
221 Efac = exp(-0.125*dt*ivxi[j+1]);
222 ivxi[j] = Efac*(ivxi[j]*Efac + 0.25*dt*GQ[j]);
225 GQ[j+1] = iQinv[j+1]*((gmx::square(ivxi[j])/iQinv[j])-kT);
232 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
239 static void boxv_trotter(t_inputrec *ir, real *veta, real dt, tensor box,
240 gmx_ekindata_t *ekind, tensor vir, real pcorr, t_extmass *MassQ)
247 tensor ekinmod, localpres;
249 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
250 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
253 if (ir->epct == epctSEMIISOTROPIC)
262 /* eta is in pure units. veta is in units of ps^-1. GW is in
263 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
264 taken to use only RATIOS of eta in updating the volume. */
266 /* we take the partial pressure tensors, modify the
267 kinetic energy tensor, and recovert to pressure */
269 if (ir->opts.nrdf[0] == 0)
271 gmx_fatal(FARGS, "Barostat is coupled to a T-group with no degrees of freedom\n");
273 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
274 alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
275 alpha *= ekind->tcstat[0].ekinscalef_nhc;
276 msmul(ekind->ekin, alpha, ekinmod);
277 /* for now, we use Elr = 0, because if you want to get it right, you
278 really should be using PME. Maybe print a warning? */
280 pscal = calc_pres(ir->ePBC, nwall, box, ekinmod, vir, localpres)+pcorr;
283 GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
289 * This file implements temperature and pressure coupling algorithms:
290 * For now only the Weak coupling and the modified weak coupling.
292 * Furthermore computation of pressure and temperature is done here
296 real calc_pres(int ePBC, int nwall, matrix box, tensor ekin, tensor vir,
302 if (ePBC == epbcNONE || (ePBC == epbcXY && nwall != 2))
308 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
309 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
313 fac = PRESFAC*2.0/det(box);
314 for (n = 0; (n < DIM); n++)
316 for (m = 0; (m < DIM); m++)
318 pres[n][m] = (ekin[n][m] - vir[n][m])*fac;
324 pr_rvecs(debug, 0, "PC: pres", pres, DIM);
325 pr_rvecs(debug, 0, "PC: ekin", ekin, DIM);
326 pr_rvecs(debug, 0, "PC: vir ", vir, DIM);
327 pr_rvecs(debug, 0, "PC: box ", box, DIM);
330 return trace(pres)/DIM;
333 real calc_temp(real ekin, real nrdf)
337 return (2.0*ekin)/(nrdf*BOLTZ);
345 void parrinellorahman_pcoupl(FILE *fplog, gmx_int64_t step,
346 const t_inputrec *ir, real dt, const tensor pres,
347 tensor box, tensor box_rel, tensor boxv,
348 tensor M, matrix mu, gmx_bool bFirstStep)
350 /* This doesn't do any coordinate updating. It just
351 * integrates the box vector equations from the calculated
352 * acceleration due to pressure difference. We also compute
353 * the tensor M which is used in update to couple the particle
354 * coordinates to the box vectors.
356 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
359 * M_nk = (h') * (h' * h + h' h) * h
361 * with the dots denoting time derivatives and h is the transformation from
362 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
363 * This also goes for the pressure and M tensors - they are transposed relative
364 * to ours. Our equation thus becomes:
367 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
369 * where b is the gromacs box matrix.
370 * Our box accelerations are given by
372 * b = vol/W inv(box') * (P-ref_P) (=h')
377 real vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
378 real atot, arel, change, maxchange, xy_pressure;
379 tensor invbox, pdiff, t1, t2;
383 gmx::invertBoxMatrix(box, invbox);
387 /* Note that PRESFAC does not occur here.
388 * The pressure and compressibility always occur as a product,
389 * therefore the pressure unit drops out.
391 maxl = std::max(box[XX][XX], box[YY][YY]);
392 maxl = std::max(maxl, box[ZZ][ZZ]);
393 for (d = 0; d < DIM; d++)
395 for (n = 0; n < DIM; n++)
398 (4*M_PI*M_PI*ir->compress[d][n])/(3*ir->tau_p*ir->tau_p*maxl);
402 m_sub(pres, ir->ref_p, pdiff);
404 if (ir->epct == epctSURFACETENSION)
406 /* Unlike Berendsen coupling it might not be trivial to include a z
407 * pressure correction here? On the other hand we don't scale the
408 * box momentarily, but change accelerations, so it might not be crucial.
410 xy_pressure = 0.5*(pres[XX][XX]+pres[YY][YY]);
411 for (d = 0; d < ZZ; d++)
413 pdiff[d][d] = (xy_pressure-(pres[ZZ][ZZ]-ir->ref_p[d][d]/box[d][d]));
417 tmmul(invbox, pdiff, t1);
418 /* Move the off-diagonal elements of the 'force' to one side to ensure
419 * that we obey the box constraints.
421 for (d = 0; d < DIM; d++)
423 for (n = 0; n < d; n++)
425 t1[d][n] += t1[n][d];
432 case epctANISOTROPIC:
433 for (d = 0; d < DIM; d++)
435 for (n = 0; n <= d; n++)
437 t1[d][n] *= winv[d][n]*vol;
442 /* calculate total volume acceleration */
443 atot = box[XX][XX]*box[YY][YY]*t1[ZZ][ZZ]+
444 box[XX][XX]*t1[YY][YY]*box[ZZ][ZZ]+
445 t1[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
447 /* set all RELATIVE box accelerations equal, and maintain total V
449 for (d = 0; d < DIM; d++)
451 for (n = 0; n <= d; n++)
453 t1[d][n] = winv[0][0]*vol*arel*box[d][n];
457 case epctSEMIISOTROPIC:
458 case epctSURFACETENSION:
459 /* Note the correction to pdiff above for surftens. coupling */
461 /* calculate total XY volume acceleration */
462 atot = box[XX][XX]*t1[YY][YY]+t1[XX][XX]*box[YY][YY];
463 arel = atot/(2*box[XX][XX]*box[YY][YY]);
464 /* set RELATIVE XY box accelerations equal, and maintain total V
465 * change speed. Dont change the third box vector accelerations */
466 for (d = 0; d < ZZ; d++)
468 for (n = 0; n <= d; n++)
470 t1[d][n] = winv[d][n]*vol*arel*box[d][n];
473 for (n = 0; n < DIM; n++)
475 t1[ZZ][n] *= winv[d][n]*vol;
479 gmx_fatal(FARGS, "Parrinello-Rahman pressure coupling type %s "
480 "not supported yet\n", EPCOUPLTYPETYPE(ir->epct));
485 for (d = 0; d < DIM; d++)
487 for (n = 0; n <= d; n++)
489 boxv[d][n] += dt*t1[d][n];
491 /* We do NOT update the box vectors themselves here, since
492 * we need them for shifting later. It is instead done last
493 * in the update() routine.
496 /* Calculate the change relative to diagonal elements-
497 since it's perfectly ok for the off-diagonal ones to
498 be zero it doesn't make sense to check the change relative
502 change = fabs(dt*boxv[d][n]/box[d][d]);
504 if (change > maxchange)
511 if (maxchange > 0.01 && fplog)
515 "\nStep %s Warning: Pressure scaling more than 1%%. "
516 "This may mean your system\n is not yet equilibrated. "
517 "Use of Parrinello-Rahman pressure coupling during\n"
518 "equilibration can lead to simulation instability, "
519 "and is discouraged.\n",
520 gmx_step_str(step, buf));
524 preserve_box_shape(ir, box_rel, boxv);
526 mtmul(boxv, box, t1); /* t1=boxv * b' */
527 mmul(invbox, t1, t2);
528 mtmul(t2, invbox, M);
530 /* Determine the scaling matrix mu for the coordinates */
531 for (d = 0; d < DIM; d++)
533 for (n = 0; n <= d; n++)
535 t1[d][n] = box[d][n] + dt*boxv[d][n];
538 preserve_box_shape(ir, box_rel, t1);
539 /* t1 is the box at t+dt, determine mu as the relative change */
540 mmul_ur0(invbox, t1, mu);
543 void berendsen_pcoupl(FILE *fplog, gmx_int64_t step,
544 const t_inputrec *ir, real dt,
545 const tensor pres, const matrix box,
549 real scalar_pressure, xy_pressure, p_corr_z;
553 * Calculate the scaling matrix mu
557 for (d = 0; d < DIM; d++)
559 scalar_pressure += pres[d][d]/DIM;
562 xy_pressure += pres[d][d]/(DIM-1);
565 /* Pressure is now in bar, everywhere. */
566 #define factor(d, m) (ir->compress[d][m]*dt/ir->tau_p)
568 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
569 * necessary for triclinic scaling
575 for (d = 0; d < DIM; d++)
577 mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
580 case epctSEMIISOTROPIC:
581 for (d = 0; d < ZZ; d++)
583 mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
586 1.0 - factor(ZZ, ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
588 case epctANISOTROPIC:
589 for (d = 0; d < DIM; d++)
591 for (n = 0; n < DIM; n++)
593 mu[d][n] = (d == n ? 1.0 : 0.0)
594 -factor(d, n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
598 case epctSURFACETENSION:
599 /* ir->ref_p[0/1] is the reference surface-tension times *
600 * the number of surfaces */
601 if (ir->compress[ZZ][ZZ])
603 p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
607 /* when the compressibity is zero, set the pressure correction *
608 * in the z-direction to zero to get the correct surface tension */
611 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
612 for (d = 0; d < DIM-1; d++)
614 mu[d][d] = 1.0 + factor(d, d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
615 - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
619 gmx_fatal(FARGS, "Berendsen pressure coupling type %s not supported yet\n",
620 EPCOUPLTYPETYPE(ir->epct));
623 /* To fullfill the orientation restrictions on triclinic boxes
624 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
625 * the other elements of mu to first order.
627 mu[YY][XX] += mu[XX][YY];
628 mu[ZZ][XX] += mu[XX][ZZ];
629 mu[ZZ][YY] += mu[YY][ZZ];
636 pr_rvecs(debug, 0, "PC: pres ", pres, 3);
637 pr_rvecs(debug, 0, "PC: mu ", mu, 3);
640 if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 ||
641 mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01 ||
642 mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
645 sprintf(buf, "\nStep %s Warning: pressure scaling more than 1%%, "
647 gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
650 fprintf(fplog, "%s", buf);
652 fprintf(stderr, "%s", buf);
656 void berendsen_pscale(const t_inputrec *ir, const matrix mu,
657 matrix box, matrix box_rel,
658 int start, int nr_atoms,
659 rvec x[], const unsigned short cFREEZE[],
662 ivec *nFreeze = ir->opts.nFreeze;
664 int nthreads gmx_unused;
666 #ifndef __clang_analyzer__
667 // cppcheck-suppress unreadVariable
668 nthreads = gmx_omp_nthreads_get(emntUpdate);
671 /* Scale the positions */
672 #pragma omp parallel for num_threads(nthreads) schedule(static)
673 for (n = start; n < start+nr_atoms; n++)
675 // Trivial OpenMP region that does not throw
689 x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
693 x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
697 x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
700 /* compute final boxlengths */
701 for (d = 0; d < DIM; d++)
703 box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
704 box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
705 box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
708 preserve_box_shape(ir, box_rel, box);
710 /* (un)shifting should NOT be done after this,
711 * since the box vectors might have changed
713 inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
716 void berendsen_tcoupl(t_inputrec *ir, gmx_ekindata_t *ekind, real dt)
720 real T, reft = 0, lll;
724 for (i = 0; (i < opts->ngtc); i++)
728 T = ekind->tcstat[i].T;
732 T = ekind->tcstat[i].Th;
735 if ((opts->tau_t[i] > 0) && (T > 0.0))
737 reft = std::max<real>(0, opts->ref_t[i]);
738 lll = std::sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
739 ekind->tcstat[i].lambda = std::max<real>(std::min<real>(lll, 1.25), 0.8);
743 ekind->tcstat[i].lambda = 1.0;
748 fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n",
749 i, T, ekind->tcstat[i].lambda);
754 void andersen_tcoupl(t_inputrec *ir, gmx_int64_t step,
755 const t_commrec *cr, const t_mdatoms *md, t_state *state, real rate, const gmx_bool *randomize, const real *boltzfac)
757 const int *gatindex = (DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
760 gmx::ThreeFry2x64<0> rng(ir->andersen_seed, gmx::RandomDomain::Thermostat);
761 gmx::UniformRealDistribution<real> uniformDist;
762 gmx::TabulatedNormalDistribution<real, 14> normalDist;
764 /* randomize the velocities of the selected particles */
766 for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */
768 int ng = gatindex ? gatindex[i] : i;
771 rng.restart(step, ng);
775 gc = md->cTC[i]; /* assign the atom to a temperature group if there are more than one */
779 if (ir->etc == etcANDERSENMASSIVE)
781 /* Randomize particle always */
786 /* Randomize particle probabilistically */
788 bRandomize = uniformDist(rng) < rate;
795 scal = std::sqrt(boltzfac[gc]*md->invmass[i]);
799 for (d = 0; d < DIM; d++)
801 state->v[i][d] = scal*normalDist(rng);
809 void nosehoover_tcoupl(t_grpopts *opts, gmx_ekindata_t *ekind, real dt,
810 double xi[], double vxi[], t_extmass *MassQ)
815 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
817 for (i = 0; (i < opts->ngtc); i++)
819 reft = std::max<real>(0, opts->ref_t[i]);
821 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
822 xi[i] += dt*(oldvxi + vxi[i])*0.5;
826 void trotter_update(t_inputrec *ir, gmx_int64_t step, gmx_ekindata_t *ekind,
827 gmx_enerdata_t *enerd, t_state *state,
828 tensor vir, t_mdatoms *md,
829 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
832 int n, i, d, ngtc, gc = 0, t;
833 t_grp_tcstat *tcstat;
835 gmx_int64_t step_eff;
837 double *scalefac, dtc;
839 rvec sumv = {0, 0, 0};
842 if (trotter_seqno <= ettTSEQ2)
844 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
845 is actually the last half step from the previous step. Thus the first half step
846 actually corresponds to the n-1 step*/
854 bCouple = (ir->nsttcouple == 1 ||
855 do_per_step(step_eff+ir->nsttcouple, ir->nsttcouple));
857 trotter_seq = trotter_seqlist[trotter_seqno];
859 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
863 dtc = ir->nsttcouple*ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
864 opts = &(ir->opts); /* just for ease of referencing */
867 snew(scalefac, opts->ngtc);
868 for (i = 0; i < ngtc; i++)
872 /* execute the series of trotter updates specified in the trotterpart array */
874 for (i = 0; i < NTROTTERPARTS; i++)
876 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
877 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
886 switch (trotter_seq[i])
890 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir,
891 enerd->term[F_PDISPCORR], MassQ);
895 NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi.data(),
896 state->nhpres_vxi.data(), NULL, &(state->veta), MassQ, FALSE);
900 NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi.data(),
901 state->nosehoover_vxi.data(), scalefac, NULL, MassQ, (ir->eI == eiVV));
902 /* need to rescale the kinetic energies and velocities here. Could
903 scale the velocities later, but we need them scaled in order to
904 produce the correct outputs, so we'll scale them here. */
906 for (t = 0; t < ngtc; t++)
908 tcstat = &ekind->tcstat[t];
909 tcstat->vscale_nhc = scalefac[t];
910 tcstat->ekinscaleh_nhc *= (scalefac[t]*scalefac[t]);
911 tcstat->ekinscalef_nhc *= (scalefac[t]*scalefac[t]);
913 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
914 /* but do we actually need the total? */
916 /* modify the velocities as well */
917 for (n = 0; n < md->homenr; n++)
919 if (md->cTC) /* does this conditional need to be here? is this always true?*/
923 for (d = 0; d < DIM; d++)
925 state->v[n][d] *= scalefac[gc];
930 for (d = 0; d < DIM; d++)
932 sumv[d] += (state->v[n][d])/md->invmass[n];
941 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
947 for (d = 0; d < DIM; d++)
949 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
951 fprintf(debug, "Conserved kappa: %15.8f %15.8f %15.8f\n", consk[0], consk[1], consk[2]);
959 extern void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit)
961 int n, i, j, d, ngtc, nh;
963 real reft, kT, ndj, nd;
965 opts = &(ir->opts); /* just for ease of referencing */
966 ngtc = ir->opts.ngtc;
967 nh = state->nhchainlength;
973 snew(MassQ->Qinv, ngtc);
975 for (i = 0; (i < ngtc); i++)
977 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
979 MassQ->Qinv[i] = 1.0/(gmx::square(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
983 MassQ->Qinv[i] = 0.0;
987 else if (EI_VV(ir->eI))
989 /* Set pressure variables */
993 if (state->vol0 == 0)
995 state->vol0 = det(state->box);
996 /* because we start by defining a fixed
997 compressibility, we need the volume at this
998 compressibility to solve the problem. */
1002 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1003 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1004 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*gmx::square(ir->tau_p/M_2PI));
1005 /* An alternate mass definition, from Tuckerman et al. */
1006 /* MassQ->Winv = 1.0/(gmx::square(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1007 for (d = 0; d < DIM; d++)
1009 for (n = 0; n < DIM; n++)
1011 MassQ->Winvm[d][n] = PRESFAC*ir->compress[d][n]/(state->vol0*gmx::square(ir->tau_p/M_2PI));
1012 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1013 before using MTTK for anisotropic states.*/
1016 /* Allocate space for thermostat variables */
1019 snew(MassQ->Qinv, ngtc*nh);
1022 /* now, set temperature variables */
1023 for (i = 0; i < ngtc; i++)
1025 if (opts->tau_t[i] > 0 && opts->ref_t[i] > 0 && opts->nrdf[i] > 0)
1027 reft = std::max<real>(0, opts->ref_t[i]);
1030 for (j = 0; j < nh; j++)
1040 MassQ->Qinv[i*nh+j] = 1.0/(gmx::square(opts->tau_t[i]/M_2PI)*ndj*kT);
1045 for (j = 0; j < nh; j++)
1047 MassQ->Qinv[i*nh+j] = 0.0;
1054 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
1056 int i, j, nnhpres, nh;
1058 real bmass, qmass, reft, kT;
1061 opts = &(ir->opts); /* just for ease of referencing */
1062 nnhpres = state->nnhpres;
1063 nh = state->nhchainlength;
1065 if (EI_VV(ir->eI) && (ir->epc == epcMTTK) && (ir->etc != etcNOSEHOOVER))
1067 gmx_fatal(FARGS, "Cannot do MTTK pressure coupling without Nose-Hoover temperature control");
1070 init_npt_masses(ir, state, MassQ, TRUE);
1072 /* first, initialize clear all the trotter calls */
1073 snew(trotter_seq, ettTSEQMAX);
1074 for (i = 0; i < ettTSEQMAX; i++)
1076 snew(trotter_seq[i], NTROTTERPARTS);
1077 for (j = 0; j < NTROTTERPARTS; j++)
1079 trotter_seq[i][j] = etrtNONE;
1081 trotter_seq[i][0] = etrtSKIPALL;
1086 /* no trotter calls, so we never use the values in the array.
1087 * We access them (so we need to define them, but ignore
1093 /* compute the kinetic energy by using the half step velocities or
1094 * the kinetic energies, depending on the order of the trotter calls */
1098 if (inputrecNptTrotter(ir))
1100 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1101 We start with the initial one. */
1102 /* first, a round that estimates veta. */
1103 trotter_seq[0][0] = etrtBAROV;
1105 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1107 /* The first half trotter update */
1108 trotter_seq[2][0] = etrtBAROV;
1109 trotter_seq[2][1] = etrtNHC;
1110 trotter_seq[2][2] = etrtBARONHC;
1112 /* The second half trotter update */
1113 trotter_seq[3][0] = etrtBARONHC;
1114 trotter_seq[3][1] = etrtNHC;
1115 trotter_seq[3][2] = etrtBAROV;
1117 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1120 else if (inputrecNvtTrotter(ir))
1122 /* This is the easy version - there are only two calls, both the same.
1123 Otherwise, even easier -- no calls */
1124 trotter_seq[2][0] = etrtNHC;
1125 trotter_seq[3][0] = etrtNHC;
1127 else if (inputrecNphTrotter(ir))
1129 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1130 We start with the initial one. */
1131 /* first, a round that estimates veta. */
1132 trotter_seq[0][0] = etrtBAROV;
1134 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1136 /* The first half trotter update */
1137 trotter_seq[2][0] = etrtBAROV;
1138 trotter_seq[2][1] = etrtBARONHC;
1140 /* The second half trotter update */
1141 trotter_seq[3][0] = etrtBARONHC;
1142 trotter_seq[3][1] = etrtBAROV;
1144 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1147 else if (ir->eI == eiVVAK)
1149 if (inputrecNptTrotter(ir))
1151 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1152 We start with the initial one. */
1153 /* first, a round that estimates veta. */
1154 trotter_seq[0][0] = etrtBAROV;
1156 /* The first half trotter update, part 1 -- double update, because it commutes */
1157 trotter_seq[1][0] = etrtNHC;
1159 /* The first half trotter update, part 2 */
1160 trotter_seq[2][0] = etrtBAROV;
1161 trotter_seq[2][1] = etrtBARONHC;
1163 /* The second half trotter update, part 1 */
1164 trotter_seq[3][0] = etrtBARONHC;
1165 trotter_seq[3][1] = etrtBAROV;
1167 /* The second half trotter update */
1168 trotter_seq[4][0] = etrtNHC;
1170 else if (inputrecNvtTrotter(ir))
1172 /* This is the easy version - there is only one call, both the same.
1173 Otherwise, even easier -- no calls */
1174 trotter_seq[1][0] = etrtNHC;
1175 trotter_seq[4][0] = etrtNHC;
1177 else if (inputrecNphTrotter(ir))
1179 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1180 We start with the initial one. */
1181 /* first, a round that estimates veta. */
1182 trotter_seq[0][0] = etrtBAROV;
1184 /* The first half trotter update, part 1 -- leave zero */
1185 trotter_seq[1][0] = etrtNHC;
1187 /* The first half trotter update, part 2 */
1188 trotter_seq[2][0] = etrtBAROV;
1189 trotter_seq[2][1] = etrtBARONHC;
1191 /* The second half trotter update, part 1 */
1192 trotter_seq[3][0] = etrtBARONHC;
1193 trotter_seq[3][1] = etrtBAROV;
1195 /* The second half trotter update -- blank for now */
1203 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
1206 snew(MassQ->QPinv, nnhpres*opts->nhchainlength);
1208 /* barostat temperature */
1209 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1211 reft = std::max<real>(0, opts->ref_t[0]);
1213 for (i = 0; i < nnhpres; i++)
1215 for (j = 0; j < nh; j++)
1225 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(gmx::square(opts->tau_t[0]/M_2PI)*qmass*kT);
1231 for (i = 0; i < nnhpres; i++)
1233 for (j = 0; j < nh; j++)
1235 MassQ->QPinv[i*nh+j] = 0.0;
1242 static real energyNoseHoover(const t_inputrec *ir, const t_state *state, const t_extmass *MassQ)
1246 int nh = state->nhchainlength;
1248 for (int i = 0; i < ir->opts.ngtc; i++)
1250 const double *ixi = &state->nosehoover_xi[i*nh];
1251 const double *ivxi = &state->nosehoover_vxi[i*nh];
1252 const double *iQinv = &(MassQ->Qinv[i*nh]);
1254 int nd = ir->opts.nrdf[i];
1255 real reft = std::max<real>(ir->opts.ref_t[i], 0);
1256 real kT = BOLTZ * reft;
1260 if (inputrecNvtTrotter(ir))
1262 /* contribution from the thermal momenta of the NH chain */
1263 for (int j = 0; j < nh; j++)
1267 energy += 0.5*gmx::square(ivxi[j])/iQinv[j];
1268 /* contribution from the thermal variable of the NH chain */
1278 energy += ndj*ixi[j]*kT;
1282 else /* Other non Trotter temperature NH control -- no chains yet. */
1284 energy += 0.5*BOLTZ*nd*gmx::square(ivxi[0])/iQinv[0];
1285 energy += nd*ixi[0]*kT;
1293 /* Returns the energy from the barostat thermostat chain */
1294 static real energyPressureMTTK(const t_inputrec *ir, const t_state *state, const t_extmass *MassQ)
1298 int nh = state->nhchainlength;
1300 for (int i = 0; i < state->nnhpres; i++)
1302 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1303 real reft = std::max<real>(ir->opts.ref_t[0], 0.0); /* using 'System' temperature */
1304 real kT = BOLTZ * reft;
1306 for (int j = 0; j < nh; j++)
1308 double iQinv = MassQ->QPinv[i*nh + j];
1311 energy += 0.5*gmx::square(state->nhpres_vxi[i*nh + j]/iQinv);
1312 /* contribution from the thermal variable of the NH chain */
1313 energy += state->nhpres_xi[i*nh + j]*kT;
1317 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]);
1325 /* Returns the energy accumulated by the V-rescale or Berendsen thermostat */
1326 static real energyVrescale(const t_inputrec *ir, const t_state *state)
1329 for (int i = 0; i < ir->opts.ngtc; i++)
1331 energy += state->therm_integral[i];
1337 real NPT_energy(const t_inputrec *ir, const t_state *state, const t_extmass *MassQ)
1341 if (ir->epc != epcNO)
1343 /* Compute the contribution of the pressure to the conserved quantity*/
1345 real vol = det(state->box);
1349 case epcPARRINELLORAHMAN:
1353 /* contribution from the pressure momenta */
1354 energyNPT += 0.5*gmx::square(state->veta)/MassQ->Winv;
1356 /* contribution from the PV term */
1357 energyNPT += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1359 if (ir->epc == epcMTTK)
1361 /* contribution from the MTTK chain */
1362 energyNPT += energyPressureMTTK(ir, state, MassQ);
1366 // Not supported, excluded in integratorHasConservedEnergyQuantity()
1370 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().");
1380 energyNPT += energyVrescale(ir, state);
1383 energyNPT += energyNoseHoover(ir, state, MassQ);
1386 case etcANDERSENMASSIVE:
1387 // Not supported, excluded in integratorHasConservedEnergyQuantity()
1390 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().");
1397 static real vrescale_sumnoises(real nn,
1398 gmx::ThreeFry2x64<> *rng,
1399 gmx::NormalDistribution<real> *normalDist)
1402 * Returns the sum of nn independent gaussian noises squared
1403 * (i.e. equivalent to summing the square of the return values
1404 * of nn calls to a normal distribution).
1406 const real ndeg_tol = 0.0001;
1408 gmx::GammaDistribution<real> gammaDist(0.5*nn, 1.0);
1410 if (nn < 2 + ndeg_tol)
1415 nn_int = (int)(nn + 0.5);
1417 if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1419 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);
1423 for (i = 0; i < nn_int; i++)
1425 gauss = (*normalDist)(*rng);
1431 /* Use a gamma distribution for any real nn > 2 */
1432 r = 2.0*gammaDist(*rng);
1438 static real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut,
1439 gmx_int64_t step, gmx_int64_t seed)
1442 * Generates a new value for the kinetic energy,
1443 * according to Bussi et al JCP (2007), Eq. (A7)
1444 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1445 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1446 * ndeg: number of degrees of freedom of the atoms to be thermalized
1447 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1449 real factor, rr, ekin_new;
1450 gmx::ThreeFry2x64<64> rng(seed, gmx::RandomDomain::Thermostat);
1451 gmx::NormalDistribution<real> normalDist;
1455 factor = exp(-1.0/taut);
1462 rng.restart(step, 0);
1464 rr = normalDist(rng);
1468 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, &rng, &normalDist) + rr*rr)/ndeg - kk) +
1469 2.0*rr*std::sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1474 void vrescale_tcoupl(t_inputrec *ir, gmx_int64_t step,
1475 gmx_ekindata_t *ekind, real dt,
1476 double therm_integral[])
1480 real Ek, Ek_ref1, Ek_ref, Ek_new;
1484 for (i = 0; (i < opts->ngtc); i++)
1488 Ek = trace(ekind->tcstat[i].ekinf);
1492 Ek = trace(ekind->tcstat[i].ekinh);
1495 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1497 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1498 Ek_ref = Ek_ref1*opts->nrdf[i];
1500 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i],
1504 /* Analytically Ek_new>=0, but we check for rounding errors */
1507 ekind->tcstat[i].lambda = 0.0;
1511 ekind->tcstat[i].lambda = std::sqrt(Ek_new/Ek);
1514 therm_integral[i] -= Ek_new - Ek;
1518 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1519 i, Ek_ref, Ek, Ek_new, ekind->tcstat[i].lambda);
1524 ekind->tcstat[i].lambda = 1.0;
1529 void rescale_velocities(gmx_ekindata_t *ekind, t_mdatoms *mdatoms,
1530 int start, int end, rvec v[])
1533 t_grp_tcstat *tcstat;
1534 unsigned short *cACC, *cTC;
1539 tcstat = ekind->tcstat;
1544 gstat = ekind->grpstat;
1545 cACC = mdatoms->cACC;
1549 for (n = start; n < end; n++)
1559 /* Only scale the velocity component relative to the COM velocity */
1560 rvec_sub(v[n], gstat[ga].u, vrel);
1561 lg = tcstat[gt].lambda;
1562 for (d = 0; d < DIM; d++)
1564 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1571 for (n = start; n < end; n++)
1577 lg = tcstat[gt].lambda;
1578 for (d = 0; d < DIM; d++)
1587 /* set target temperatures if we are annealing */
1588 void update_annealing_target_temp(t_inputrec *ir, real t, gmx_update_t *upd)
1590 int i, j, n, npoints;
1591 real pert, thist = 0, x;
1593 for (i = 0; i < ir->opts.ngtc; i++)
1595 npoints = ir->opts.anneal_npoints[i];
1596 switch (ir->opts.annealing[i])
1601 /* calculate time modulo the period */
1602 pert = ir->opts.anneal_time[i][npoints-1];
1603 n = static_cast<int>(t / pert);
1604 thist = t - n*pert; /* modulo time */
1605 /* Make sure rounding didn't get us outside the interval */
1606 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1615 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)", i, ir->opts.ngtc, npoints);
1617 /* We are doing annealing for this group if we got here,
1618 * and we have the (relative) time as thist.
1619 * calculate target temp */
1621 while ((j < npoints-1) && (thist > (ir->opts.anneal_time[i][j+1])))
1627 /* Found our position between points j and j+1.
1628 * Interpolate: x is the amount from j+1, (1-x) from point j
1629 * First treat possible jumps in temperature as a special case.
1631 if ((ir->opts.anneal_time[i][j+1]-ir->opts.anneal_time[i][j]) < GMX_REAL_EPS*100)
1633 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][j+1];
1637 x = ((thist-ir->opts.anneal_time[i][j])/
1638 (ir->opts.anneal_time[i][j+1]-ir->opts.anneal_time[i][j]));
1639 ir->opts.ref_t[i] = x*ir->opts.anneal_temp[i][j+1]+(1-x)*ir->opts.anneal_temp[i][j];
1644 ir->opts.ref_t[i] = ir->opts.anneal_temp[i][npoints-1];
1648 update_temperature_constants(upd, ir);