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, 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.
49 #include "gmx_fatal.h"
52 #include "gromacs/random/random.h"
56 #define NTROTTERPARTS 3
58 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
60 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
61 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
63 #define MAX_SUZUKI_YOSHIDA_NUM 5
64 #define SUZUKI_YOSHIDA_NUM 5
66 static const double sy_const_1[] = { 1. };
67 static const double sy_const_3[] = { 0.828981543588751, -0.657963087177502, 0.828981543588751 };
68 static const double sy_const_5[] = { 0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065, 0.2967324292201065 };
70 static const double* sy_const[] = {
80 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
84 {0.828981543588751,-0.657963087177502,0.828981543588751},
86 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
89 /* these integration routines are only referenced inside this file */
90 static void NHC_trotter(t_grpopts *opts, int nvar, gmx_ekindata_t *ekind, real dtfull,
91 double xi[], double vxi[], double scalefac[], real *veta, t_extmass *MassQ, gmx_bool bEkinAveVel)
94 /* general routine for both barostat and thermostat nose hoover chains */
96 int i, j, mi, mj, jmax;
97 double Ekin, Efac, reft, kT, nd;
104 int mstepsi, mstepsj;
105 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
106 int nh = opts->nhchainlength;
109 mstepsi = mstepsj = ns;
111 /* if scalefac is NULL, we are doing the NHC of the barostat */
114 if (scalefac == NULL)
119 for (i = 0; i < nvar; i++)
122 /* make it easier to iterate by selecting
123 out the sub-array that corresponds to this T group */
129 iQinv = &(MassQ->QPinv[i*nh]);
130 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
131 reft = max(0.0, opts->ref_t[0]);
132 Ekin = sqr(*veta)/MassQ->Winv;
136 iQinv = &(MassQ->Qinv[i*nh]);
137 tcstat = &ekind->tcstat[i];
139 reft = max(0.0, opts->ref_t[i]);
142 Ekin = 2*trace(tcstat->ekinf)*tcstat->ekinscalef_nhc;
146 Ekin = 2*trace(tcstat->ekinh)*tcstat->ekinscaleh_nhc;
151 for (mi = 0; mi < mstepsi; mi++)
153 for (mj = 0; mj < mstepsj; mj++)
155 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
156 dt = sy_const[ns][mj] * dtfull / mstepsi;
158 /* compute the thermal forces */
159 GQ[0] = iQinv[0]*(Ekin - nd*kT);
161 for (j = 0; j < nh-1; j++)
165 /* we actually don't need to update here if we save the
166 state of the GQ, but it's easier to just recompute*/
167 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
175 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
176 for (j = nh-1; j > 0; j--)
178 Efac = exp(-0.125*dt*ivxi[j]);
179 ivxi[j-1] = Efac*(ivxi[j-1]*Efac + 0.25*dt*GQ[j-1]);
182 Efac = exp(-0.5*dt*ivxi[0]);
193 /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
194 able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
195 think about this a bit more . . . */
197 GQ[0] = iQinv[0]*(Ekin - nd*kT);
199 /* update thermostat positions */
200 for (j = 0; j < nh; j++)
202 ixi[j] += 0.5*dt*ivxi[j];
205 for (j = 0; j < nh-1; j++)
207 Efac = exp(-0.125*dt*ivxi[j+1]);
208 ivxi[j] = Efac*(ivxi[j]*Efac + 0.25*dt*GQ[j]);
211 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
218 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
225 static void boxv_trotter(t_inputrec *ir, real *veta, real dt, tensor box,
226 gmx_ekindata_t *ekind, tensor vir, real pcorr, t_extmass *MassQ)
231 int i, j, d, n, nwall;
233 tensor Winvm, ekinmod, localpres;
235 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
236 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
239 if (ir->epct == epctSEMIISOTROPIC)
248 /* eta is in pure units. veta is in units of ps^-1. GW is in
249 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
250 taken to use only RATIOS of eta in updating the volume. */
252 /* we take the partial pressure tensors, modify the
253 kinetic energy tensor, and recovert to pressure */
255 if (ir->opts.nrdf[0] == 0)
257 gmx_fatal(FARGS, "Barostat is coupled to a T-group with no degrees of freedom\n");
259 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
260 alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
261 alpha *= ekind->tcstat[0].ekinscalef_nhc;
262 msmul(ekind->ekin, alpha, ekinmod);
263 /* for now, we use Elr = 0, because if you want to get it right, you
264 really should be using PME. Maybe print a warning? */
266 pscal = calc_pres(ir->ePBC, nwall, box, ekinmod, vir, localpres)+pcorr;
269 GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
275 * This file implements temperature and pressure coupling algorithms:
276 * For now only the Weak coupling and the modified weak coupling.
278 * Furthermore computation of pressure and temperature is done here
282 real calc_pres(int ePBC, int nwall, matrix box, tensor ekin, tensor vir,
288 if (ePBC == epbcNONE || (ePBC == epbcXY && nwall != 2))
294 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
295 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
299 fac = PRESFAC*2.0/det(box);
300 for (n = 0; (n < DIM); n++)
302 for (m = 0; (m < DIM); m++)
304 pres[n][m] = (ekin[n][m] - vir[n][m])*fac;
310 pr_rvecs(debug, 0, "PC: pres", pres, DIM);
311 pr_rvecs(debug, 0, "PC: ekin", ekin, DIM);
312 pr_rvecs(debug, 0, "PC: vir ", vir, DIM);
313 pr_rvecs(debug, 0, "PC: box ", box, DIM);
316 return trace(pres)/DIM;
319 real calc_temp(real ekin, real nrdf)
323 return (2.0*ekin)/(nrdf*BOLTZ);
331 void parrinellorahman_pcoupl(FILE *fplog, gmx_int64_t step,
332 t_inputrec *ir, real dt, tensor pres,
333 tensor box, tensor box_rel, tensor boxv,
334 tensor M, matrix mu, gmx_bool bFirstStep)
336 /* This doesn't do any coordinate updating. It just
337 * integrates the box vector equations from the calculated
338 * acceleration due to pressure difference. We also compute
339 * the tensor M which is used in update to couple the particle
340 * coordinates to the box vectors.
342 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
345 * M_nk = (h') * (h' * h + h' h) * h
347 * with the dots denoting time derivatives and h is the transformation from
348 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
349 * This also goes for the pressure and M tensors - they are transposed relative
350 * to ours. Our equation thus becomes:
353 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
355 * where b is the gromacs box matrix.
356 * Our box accelerations are given by
358 * b = vol/W inv(box') * (P-ref_P) (=h')
363 real vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
364 real atot, arel, change, maxchange, xy_pressure;
365 tensor invbox, pdiff, t1, t2;
369 m_inv_ur0(box, invbox);
373 /* Note that PRESFAC does not occur here.
374 * The pressure and compressibility always occur as a product,
375 * therefore the pressure unit drops out.
377 maxl = max(box[XX][XX], box[YY][YY]);
378 maxl = max(maxl, box[ZZ][ZZ]);
379 for (d = 0; d < DIM; d++)
381 for (n = 0; n < DIM; n++)
384 (4*M_PI*M_PI*ir->compress[d][n])/(3*ir->tau_p*ir->tau_p*maxl);
388 m_sub(pres, ir->ref_p, pdiff);
390 if (ir->epct == epctSURFACETENSION)
392 /* Unlike Berendsen coupling it might not be trivial to include a z
393 * pressure correction here? On the other hand we don't scale the
394 * box momentarily, but change accelerations, so it might not be crucial.
396 xy_pressure = 0.5*(pres[XX][XX]+pres[YY][YY]);
397 for (d = 0; d < ZZ; d++)
399 pdiff[d][d] = (xy_pressure-(pres[ZZ][ZZ]-ir->ref_p[d][d]/box[d][d]));
403 tmmul(invbox, pdiff, t1);
404 /* Move the off-diagonal elements of the 'force' to one side to ensure
405 * that we obey the box constraints.
407 for (d = 0; d < DIM; d++)
409 for (n = 0; n < d; n++)
411 t1[d][n] += t1[n][d];
418 case epctANISOTROPIC:
419 for (d = 0; d < DIM; d++)
421 for (n = 0; n <= d; n++)
423 t1[d][n] *= winv[d][n]*vol;
428 /* calculate total volume acceleration */
429 atot = box[XX][XX]*box[YY][YY]*t1[ZZ][ZZ]+
430 box[XX][XX]*t1[YY][YY]*box[ZZ][ZZ]+
431 t1[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
433 /* set all RELATIVE box accelerations equal, and maintain total V
435 for (d = 0; d < DIM; d++)
437 for (n = 0; n <= d; n++)
439 t1[d][n] = winv[0][0]*vol*arel*box[d][n];
443 case epctSEMIISOTROPIC:
444 case epctSURFACETENSION:
445 /* Note the correction to pdiff above for surftens. coupling */
447 /* calculate total XY volume acceleration */
448 atot = box[XX][XX]*t1[YY][YY]+t1[XX][XX]*box[YY][YY];
449 arel = atot/(2*box[XX][XX]*box[YY][YY]);
450 /* set RELATIVE XY box accelerations equal, and maintain total V
451 * change speed. Dont change the third box vector accelerations */
452 for (d = 0; d < ZZ; d++)
454 for (n = 0; n <= d; n++)
456 t1[d][n] = winv[d][n]*vol*arel*box[d][n];
459 for (n = 0; n < DIM; n++)
461 t1[ZZ][n] *= winv[d][n]*vol;
465 gmx_fatal(FARGS, "Parrinello-Rahman pressure coupling type %s "
466 "not supported yet\n", EPCOUPLTYPETYPE(ir->epct));
471 for (d = 0; d < DIM; d++)
473 for (n = 0; n <= d; n++)
475 boxv[d][n] += dt*t1[d][n];
477 /* We do NOT update the box vectors themselves here, since
478 * we need them for shifting later. It is instead done last
479 * in the update() routine.
482 /* Calculate the change relative to diagonal elements-
483 since it's perfectly ok for the off-diagonal ones to
484 be zero it doesn't make sense to check the change relative
488 change = fabs(dt*boxv[d][n]/box[d][d]);
490 if (change > maxchange)
497 if (maxchange > 0.01 && fplog)
501 "\nStep %s Warning: Pressure scaling more than 1%%. "
502 "This may mean your system\n is not yet equilibrated. "
503 "Use of Parrinello-Rahman pressure coupling during\n"
504 "equilibration can lead to simulation instability, "
505 "and is discouraged.\n",
506 gmx_step_str(step, buf));
510 preserve_box_shape(ir, box_rel, boxv);
512 mtmul(boxv, box, t1); /* t1=boxv * b' */
513 mmul(invbox, t1, t2);
514 mtmul(t2, invbox, M);
516 /* Determine the scaling matrix mu for the coordinates */
517 for (d = 0; d < DIM; d++)
519 for (n = 0; n <= d; n++)
521 t1[d][n] = box[d][n] + dt*boxv[d][n];
524 preserve_box_shape(ir, box_rel, t1);
525 /* t1 is the box at t+dt, determine mu as the relative change */
526 mmul_ur0(invbox, t1, mu);
529 void berendsen_pcoupl(FILE *fplog, gmx_int64_t step,
530 t_inputrec *ir, real dt, tensor pres, matrix box,
534 real scalar_pressure, xy_pressure, p_corr_z;
535 char *ptr, buf[STRLEN];
538 * Calculate the scaling matrix mu
542 for (d = 0; d < DIM; d++)
544 scalar_pressure += pres[d][d]/DIM;
547 xy_pressure += pres[d][d]/(DIM-1);
550 /* Pressure is now in bar, everywhere. */
551 #define factor(d, m) (ir->compress[d][m]*dt/ir->tau_p)
553 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
554 * necessary for triclinic scaling
560 for (d = 0; d < DIM; d++)
562 mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
565 case epctSEMIISOTROPIC:
566 for (d = 0; d < ZZ; d++)
568 mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
571 1.0 - factor(ZZ, ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
573 case epctANISOTROPIC:
574 for (d = 0; d < DIM; d++)
576 for (n = 0; n < DIM; n++)
578 mu[d][n] = (d == n ? 1.0 : 0.0)
579 -factor(d, n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
583 case epctSURFACETENSION:
584 /* ir->ref_p[0/1] is the reference surface-tension times *
585 * the number of surfaces */
586 if (ir->compress[ZZ][ZZ])
588 p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
592 /* when the compressibity is zero, set the pressure correction *
593 * in the z-direction to zero to get the correct surface tension */
596 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
597 for (d = 0; d < DIM-1; d++)
599 mu[d][d] = 1.0 + factor(d, d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
600 - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
604 gmx_fatal(FARGS, "Berendsen pressure coupling type %s not supported yet\n",
605 EPCOUPLTYPETYPE(ir->epct));
608 /* To fullfill the orientation restrictions on triclinic boxes
609 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
610 * the other elements of mu to first order.
612 mu[YY][XX] += mu[XX][YY];
613 mu[ZZ][XX] += mu[XX][ZZ];
614 mu[ZZ][YY] += mu[YY][ZZ];
621 pr_rvecs(debug, 0, "PC: pres ", pres, 3);
622 pr_rvecs(debug, 0, "PC: mu ", mu, 3);
625 if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 ||
626 mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01 ||
627 mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
630 sprintf(buf, "\nStep %s Warning: pressure scaling more than 1%%, "
632 gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
635 fprintf(fplog, "%s", buf);
637 fprintf(stderr, "%s", buf);
641 void berendsen_pscale(t_inputrec *ir, matrix mu,
642 matrix box, matrix box_rel,
643 int start, int nr_atoms,
644 rvec x[], unsigned short cFREEZE[],
647 ivec *nFreeze = ir->opts.nFreeze;
650 /* Scale the positions */
651 for (n = start; n < start+nr_atoms; n++)
660 x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
664 x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
668 x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
671 /* compute final boxlengths */
672 for (d = 0; d < DIM; d++)
674 box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
675 box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
676 box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
679 preserve_box_shape(ir, box_rel, box);
681 /* (un)shifting should NOT be done after this,
682 * since the box vectors might have changed
684 inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
687 void berendsen_tcoupl(t_inputrec *ir, gmx_ekindata_t *ekind, real dt)
691 real T, reft = 0, lll;
695 for (i = 0; (i < opts->ngtc); i++)
699 T = ekind->tcstat[i].T;
703 T = ekind->tcstat[i].Th;
706 if ((opts->tau_t[i] > 0) && (T > 0.0))
708 reft = max(0.0, opts->ref_t[i]);
709 lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
710 ekind->tcstat[i].lambda = max(min(lll, 1.25), 0.8);
714 ekind->tcstat[i].lambda = 1.0;
719 fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n",
720 i, T, ekind->tcstat[i].lambda);
725 static int poisson_variate(real lambda, gmx_rng_t rng)
737 p *= gmx_rng_uniform_real(rng);
744 void andersen_tcoupl(t_inputrec *ir, t_mdatoms *md, t_state *state, gmx_rng_t rng, real rate, t_idef *idef, int nblocks, int *sblock, gmx_bool *randatom, int *randatom_list, gmx_bool *randomize, real *boltzfac)
747 int i, j, k, d, len, n, ngtc, gc = 0;
748 int nshake, nsettle, nrandom, nrand_group;
749 real boltz, scal, reft, prand;
752 /* convenience variables */
756 /* idef is only passed in if it's chance-per-particle andersen, so
757 it essentially serves as a boolean to determine which type of
758 andersen is being used */
762 /* randomly atoms to randomize. However, all constraint
763 groups have to have either all of the atoms or none of the
767 1. Select whether or not to randomize each atom to get the correct probability.
768 2. Cycle through the constraint groups.
769 2a. for each constraint group, determine the fraction f of that constraint group that are
770 chosen to be randomized.
771 2b. all atoms in the constraint group are randomized with probability f.
775 if ((rate < 0.05) && (md->homenr > 50))
777 /* if the rate is relatively high, use a standard method, if low rate,
779 /* poisson distributions approxmation, more efficient for
780 * low rates, fewer random numbers required */
781 nrandom = poisson_variate(md->homenr*rate, rng); /* how many do we randomize? Use Poisson. */
782 /* now we know how many, choose them randomly. No worries about repeats, at this rate, it's negligible.
783 worst thing that happens, it lowers the true rate an negligible amount */
784 for (i = 0; i < nrandom; i++)
786 randatom[(int)(gmx_rng_uniform_real(rng)*md->homenr)] = TRUE;
791 for (i = 0; i < md->homenr; i++)
793 if (gmx_rng_uniform_real(rng) < rate)
801 /* instead of looping over the constraint groups, if we had a
802 list of which atoms were in which constraint groups, we
803 could then loop over only the groups that are randomized
804 now. But that is not available now. Create later after
805 determining whether there actually is any slowing. */
807 /* first, loop through the settles to make sure all groups either entirely randomized, or not randomized. */
809 nsettle = idef->il[F_SETTLE].nr/2;
810 for (i = 0; i < nsettle; i++)
812 iatoms = idef->il[F_SETTLE].iatoms;
814 for (k = 0; k < 3; k++) /* settles are always 3 atoms, hardcoded */
816 if (randatom[iatoms[2*i+1]+k])
818 nrand_group++; /* count the number of atoms to be shaken in the settles group */
819 randatom[iatoms[2*i+1]+k] = FALSE;
825 prand = (nrand_group)/3.0; /* use this fraction to compute the probability the
826 whole group is randomized */
827 if (gmx_rng_uniform_real(rng) < prand)
829 for (k = 0; k < 3; k++)
831 randatom[iatoms[2*i+1]+k] = TRUE; /* mark them all to be randomized */
838 /* now loop through the shake groups */
840 for (i = 0; i < nshake; i++)
842 iatoms = &(idef->il[F_CONSTR].iatoms[sblock[i]]);
843 len = sblock[i+1]-sblock[i];
845 for (k = 0; k < len; k++)
848 { /* only 2/3 of the sblock items are atoms, the others are labels */
849 if (randatom[iatoms[k]])
852 randatom[iatoms[k]] = FALSE; /* need to mark it false here in case the atom is in more than
853 one group in the shake block */
860 prand = (nrand_group)/(1.0*(2*len/3));
861 if (gmx_rng_uniform_real(rng) < prand)
863 for (k = 0; k < len; k++)
866 { /* only 2/3 of the sblock items are atoms, the others are labels */
867 randatom[iatoms[k]] = TRUE; /* randomize all of them */
877 for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */
881 randatom_list[n] = i;
885 nrandom = n; /* there are some values of nrandom for which
886 this algorithm won't work; for example all
887 water molecules and nrandom =/= 3. Better to
888 recount and use this number (which we
889 calculate anyway: it will not affect
890 the average number of atoms accepted.
896 /* if it's andersen-massive, then randomize all the atoms */
897 nrandom = md->homenr;
898 for (i = 0; i < nrandom; i++)
900 randatom_list[i] = i;
904 /* randomize the velocities of the selected particles */
906 for (i = 0; i < nrandom; i++) /* now loop over the list of atoms */
908 n = randatom_list[i];
911 gc = md->cTC[n]; /* assign the atom to a temperature group if there are more than one */
915 scal = sqrt(boltzfac[gc]*md->invmass[n]);
916 for (d = 0; d < DIM; d++)
918 state->v[n][d] = scal*gmx_rng_gaussian_table(rng);
921 randatom[n] = FALSE; /* unmark this atom for randomization */
926 void nosehoover_tcoupl(t_grpopts *opts, gmx_ekindata_t *ekind, real dt,
927 double xi[], double vxi[], t_extmass *MassQ)
932 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
934 for (i = 0; (i < opts->ngtc); i++)
936 reft = max(0.0, opts->ref_t[i]);
938 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
939 xi[i] += dt*(oldvxi + vxi[i])*0.5;
943 t_state *init_bufstate(const t_state *template_state)
946 int nc = template_state->nhchainlength;
948 snew(state->nosehoover_xi, nc*template_state->ngtc);
949 snew(state->nosehoover_vxi, nc*template_state->ngtc);
950 snew(state->therm_integral, template_state->ngtc);
951 snew(state->nhpres_xi, nc*template_state->nnhpres);
952 snew(state->nhpres_vxi, nc*template_state->nnhpres);
957 void destroy_bufstate(t_state *state)
961 sfree(state->nosehoover_xi);
962 sfree(state->nosehoover_vxi);
963 sfree(state->therm_integral);
964 sfree(state->nhpres_xi);
965 sfree(state->nhpres_vxi);
969 void trotter_update(t_inputrec *ir, gmx_int64_t step, gmx_ekindata_t *ekind,
970 gmx_enerdata_t *enerd, t_state *state,
971 tensor vir, t_mdatoms *md,
972 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
975 int n, i, j, d, ntgrp, ngtc, gc = 0;
976 t_grp_tcstat *tcstat;
978 gmx_int64_t step_eff;
979 real ecorr, pcorr, dvdlcorr;
980 real bmass, qmass, reft, kT, dt, nd;
981 tensor dumpres, dumvir;
982 double *scalefac, dtc;
984 rvec sumv = {0, 0, 0}, consk;
987 if (trotter_seqno <= ettTSEQ2)
989 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
990 is actually the last half step from the previous step. Thus the first half step
991 actually corresponds to the n-1 step*/
999 bCouple = (ir->nsttcouple == 1 ||
1000 do_per_step(step_eff+ir->nsttcouple, ir->nsttcouple));
1002 trotter_seq = trotter_seqlist[trotter_seqno];
1004 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
1008 dtc = ir->nsttcouple*ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
1009 opts = &(ir->opts); /* just for ease of referencing */
1012 snew(scalefac, opts->ngtc);
1013 for (i = 0; i < ngtc; i++)
1017 /* execute the series of trotter updates specified in the trotterpart array */
1019 for (i = 0; i < NTROTTERPARTS; i++)
1021 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
1022 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
1031 switch (trotter_seq[i])
1035 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir,
1036 enerd->term[F_PDISPCORR], MassQ);
1040 NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi,
1041 state->nhpres_vxi, NULL, &(state->veta), MassQ, FALSE);
1045 NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi,
1046 state->nosehoover_vxi, scalefac, NULL, MassQ, (ir->eI == eiVV));
1047 /* need to rescale the kinetic energies and velocities here. Could
1048 scale the velocities later, but we need them scaled in order to
1049 produce the correct outputs, so we'll scale them here. */
1051 for (i = 0; i < ngtc; i++)
1053 tcstat = &ekind->tcstat[i];
1054 tcstat->vscale_nhc = scalefac[i];
1055 tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
1056 tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
1058 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
1059 /* but do we actually need the total? */
1061 /* modify the velocities as well */
1062 for (n = md->start; n < md->start+md->homenr; n++)
1064 if (md->cTC) /* does this conditional need to be here? is this always true?*/
1068 for (d = 0; d < DIM; d++)
1070 state->v[n][d] *= scalefac[gc];
1075 for (d = 0; d < DIM; d++)
1077 sumv[d] += (state->v[n][d])/md->invmass[n];
1086 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
1092 for (d = 0; d < DIM; d++)
1094 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
1096 fprintf(debug, "Conserved kappa: %15.8f %15.8f %15.8f\n", consk[0], consk[1], consk[2]);
1104 extern void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit)
1106 int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
1107 t_grp_tcstat *tcstat;
1109 real ecorr, pcorr, dvdlcorr;
1110 real bmass, qmass, reft, kT, dt, ndj, nd;
1111 tensor dumpres, dumvir;
1113 opts = &(ir->opts); /* just for ease of referencing */
1114 ngtc = ir->opts.ngtc;
1115 nnhpres = state->nnhpres;
1116 nh = state->nhchainlength;
1122 snew(MassQ->Qinv, ngtc);
1124 for (i = 0; (i < ngtc); i++)
1126 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1128 MassQ->Qinv[i] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
1132 MassQ->Qinv[i] = 0.0;
1136 else if (EI_VV(ir->eI))
1138 /* Set pressure variables */
1142 if (state->vol0 == 0)
1144 state->vol0 = det(state->box);
1145 /* because we start by defining a fixed
1146 compressibility, we need the volume at this
1147 compressibility to solve the problem. */
1151 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1152 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1153 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
1154 /* An alternate mass definition, from Tuckerman et al. */
1155 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1156 for (d = 0; d < DIM; d++)
1158 for (n = 0; n < DIM; n++)
1160 MassQ->Winvm[d][n] = PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
1161 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1162 before using MTTK for anisotropic states.*/
1165 /* Allocate space for thermostat variables */
1168 snew(MassQ->Qinv, ngtc*nh);
1171 /* now, set temperature variables */
1172 for (i = 0; i < ngtc; i++)
1174 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1176 reft = max(0.0, opts->ref_t[i]);
1179 for (j = 0; j < nh; j++)
1189 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
1195 for (j = 0; j < nh; j++)
1197 MassQ->Qinv[i*nh+j] = 0.0;
1204 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
1206 int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
1207 t_grp_tcstat *tcstat;
1209 real ecorr, pcorr, dvdlcorr;
1210 real bmass, qmass, reft, kT, dt, ndj, nd;
1211 tensor dumpres, dumvir;
1214 opts = &(ir->opts); /* just for ease of referencing */
1216 nnhpres = state->nnhpres;
1217 nh = state->nhchainlength;
1219 init_npt_masses(ir, state, MassQ, TRUE);
1221 /* first, initialize clear all the trotter calls */
1222 snew(trotter_seq, ettTSEQMAX);
1223 for (i = 0; i < ettTSEQMAX; i++)
1225 snew(trotter_seq[i], NTROTTERPARTS);
1226 for (j = 0; j < NTROTTERPARTS; j++)
1228 trotter_seq[i][j] = etrtNONE;
1230 trotter_seq[i][0] = etrtSKIPALL;
1235 /* no trotter calls, so we never use the values in the array.
1236 * We access them (so we need to define them, but ignore
1242 /* compute the kinetic energy by using the half step velocities or
1243 * the kinetic energies, depending on the order of the trotter calls */
1247 if (IR_NPT_TROTTER(ir))
1249 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1250 We start with the initial one. */
1251 /* first, a round that estimates veta. */
1252 trotter_seq[0][0] = etrtBAROV;
1254 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1256 /* The first half trotter update */
1257 trotter_seq[2][0] = etrtBAROV;
1258 trotter_seq[2][1] = etrtNHC;
1259 trotter_seq[2][2] = etrtBARONHC;
1261 /* The second half trotter update */
1262 trotter_seq[3][0] = etrtBARONHC;
1263 trotter_seq[3][1] = etrtNHC;
1264 trotter_seq[3][2] = etrtBAROV;
1266 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1269 else if (IR_NVT_TROTTER(ir))
1271 /* This is the easy version - there are only two calls, both the same.
1272 Otherwise, even easier -- no calls */
1273 trotter_seq[2][0] = etrtNHC;
1274 trotter_seq[3][0] = etrtNHC;
1276 else if (IR_NPH_TROTTER(ir))
1278 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1279 We start with the initial one. */
1280 /* first, a round that estimates veta. */
1281 trotter_seq[0][0] = etrtBAROV;
1283 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1285 /* The first half trotter update */
1286 trotter_seq[2][0] = etrtBAROV;
1287 trotter_seq[2][1] = etrtBARONHC;
1289 /* The second half trotter update */
1290 trotter_seq[3][0] = etrtBARONHC;
1291 trotter_seq[3][1] = etrtBAROV;
1293 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1296 else if (ir->eI == eiVVAK)
1298 if (IR_NPT_TROTTER(ir))
1300 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1301 We start with the initial one. */
1302 /* first, a round that estimates veta. */
1303 trotter_seq[0][0] = etrtBAROV;
1305 /* The first half trotter update, part 1 -- double update, because it commutes */
1306 trotter_seq[1][0] = etrtNHC;
1308 /* The first half trotter update, part 2 */
1309 trotter_seq[2][0] = etrtBAROV;
1310 trotter_seq[2][1] = etrtBARONHC;
1312 /* The second half trotter update, part 1 */
1313 trotter_seq[3][0] = etrtBARONHC;
1314 trotter_seq[3][1] = etrtBAROV;
1316 /* The second half trotter update */
1317 trotter_seq[4][0] = etrtNHC;
1319 else if (IR_NVT_TROTTER(ir))
1321 /* This is the easy version - there is only one call, both the same.
1322 Otherwise, even easier -- no calls */
1323 trotter_seq[1][0] = etrtNHC;
1324 trotter_seq[4][0] = etrtNHC;
1326 else if (IR_NPH_TROTTER(ir))
1328 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1329 We start with the initial one. */
1330 /* first, a round that estimates veta. */
1331 trotter_seq[0][0] = etrtBAROV;
1333 /* The first half trotter update, part 1 -- leave zero */
1334 trotter_seq[1][0] = etrtNHC;
1336 /* The first half trotter update, part 2 */
1337 trotter_seq[2][0] = etrtBAROV;
1338 trotter_seq[2][1] = etrtBARONHC;
1340 /* The second half trotter update, part 1 */
1341 trotter_seq[3][0] = etrtBARONHC;
1342 trotter_seq[3][1] = etrtBAROV;
1344 /* The second half trotter update -- blank for now */
1352 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
1355 snew(MassQ->QPinv, nnhpres*opts->nhchainlength);
1357 /* barostat temperature */
1358 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1360 reft = max(0.0, opts->ref_t[0]);
1362 for (i = 0; i < nnhpres; i++)
1364 for (j = 0; j < nh; j++)
1374 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1380 for (i = 0; i < nnhpres; i++)
1382 for (j = 0; j < nh; j++)
1384 MassQ->QPinv[i*nh+j] = 0.0;
1391 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1393 int i, j, nd, ndj, bmass, qmass, ngtcall;
1394 real ener_npt, reft, eta, kT, tau;
1397 real vol, dbaro, W, Q;
1398 int nh = state->nhchainlength;
1402 /* now we compute the contribution of the pressure to the conserved quantity*/
1404 if (ir->epc == epcMTTK)
1406 /* find the volume, and the kinetic energy of the volume */
1412 /* contribution from the pressure momenenta */
1413 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1415 /* contribution from the PV term */
1416 vol = det(state->box);
1417 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1420 case epctANISOTROPIC:
1424 case epctSURFACETENSION:
1427 case epctSEMIISOTROPIC:
1435 if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1437 /* add the energy from the barostat thermostat chain */
1438 for (i = 0; i < state->nnhpres; i++)
1441 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1442 ivxi = &state->nhpres_vxi[i*nh];
1443 ixi = &state->nhpres_xi[i*nh];
1444 iQinv = &(MassQ->QPinv[i*nh]);
1445 reft = max(ir->opts.ref_t[0], 0); /* using 'System' temperature */
1448 for (j = 0; j < nh; j++)
1452 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1453 /* contribution from the thermal variable of the NH chain */
1454 ener_npt += ixi[j]*kT;
1458 fprintf(debug, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i, j, ivxi[j], ixi[j]);
1466 for (i = 0; i < ir->opts.ngtc; i++)
1468 ixi = &state->nosehoover_xi[i*nh];
1469 ivxi = &state->nosehoover_vxi[i*nh];
1470 iQinv = &(MassQ->Qinv[i*nh]);
1472 nd = ir->opts.nrdf[i];
1473 reft = max(ir->opts.ref_t[i], 0);
1478 if (IR_NVT_TROTTER(ir))
1480 /* contribution from the thermal momenta of the NH chain */
1481 for (j = 0; j < nh; j++)
1485 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1486 /* contribution from the thermal variable of the NH chain */
1495 ener_npt += ndj*ixi[j]*kT;
1499 else /* Other non Trotter temperature NH control -- no chains yet. */
1501 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1502 ener_npt += nd*ixi[0]*kT;
1510 static real vrescale_gamdev(int ia, gmx_rng_t rng)
1511 /* Gamma distribution, adapted from numerical recipes */
1514 real am, e, s, v1, v2, x, y;
1521 for (j = 1; j <= ia; j++)
1523 x *= gmx_rng_uniform_real(rng);
1537 v1 = gmx_rng_uniform_real(rng);
1538 v2 = 2.0*gmx_rng_uniform_real(rng)-1.0;
1540 while (v1*v1 + v2*v2 > 1.0 ||
1541 v1*v1*GMX_REAL_MAX < 3.0*ia);
1542 /* The last check above ensures that both x (3.0 > 2.0 in s)
1543 * and the pre-factor for e do not go out of range.
1547 s = sqrt(2.0*am + 1.0);
1551 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1553 while (gmx_rng_uniform_real(rng) > e);
1559 static real vrescale_sumnoises(int nn, gmx_rng_t rng)
1562 * Returns the sum of n independent gaussian noises squared
1563 * (i.e. equivalent to summing the square of the return values
1564 * of nn calls to gmx_rng_gaussian_real).xs
1574 rr = gmx_rng_gaussian_real(rng);
1577 else if (nn % 2 == 0)
1579 return 2.0*vrescale_gamdev(nn/2, rng);
1583 rr = gmx_rng_gaussian_real(rng);
1584 return 2.0*vrescale_gamdev((nn-1)/2, rng) + rr*rr;
1588 static real vrescale_resamplekin(real kk, real sigma, int ndeg, real taut,
1592 * Generates a new value for the kinetic energy,
1593 * according to Bussi et al JCP (2007), Eq. (A7)
1594 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1595 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1596 * ndeg: number of degrees of freedom of the atoms to be thermalized
1597 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1603 factor = exp(-1.0/taut);
1609 rr = gmx_rng_gaussian_real(rng);
1612 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, rng) + rr*rr)/ndeg - kk) +
1613 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1616 void vrescale_tcoupl(t_inputrec *ir, gmx_ekindata_t *ekind, real dt,
1617 double therm_integral[], gmx_rng_t rng)
1621 real Ek, Ek_ref1, Ek_ref, Ek_new;
1625 for (i = 0; (i < opts->ngtc); i++)
1629 Ek = trace(ekind->tcstat[i].ekinf);
1633 Ek = trace(ekind->tcstat[i].ekinh);
1636 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1638 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1639 Ek_ref = Ek_ref1*opts->nrdf[i];
1641 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i],
1642 opts->tau_t[i]/dt, rng);
1644 /* Analytically Ek_new>=0, but we check for rounding errors */
1647 ekind->tcstat[i].lambda = 0.0;
1651 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1654 therm_integral[i] -= Ek_new - Ek;
1658 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1659 i, Ek_ref, Ek, Ek_new, ekind->tcstat[i].lambda);
1664 ekind->tcstat[i].lambda = 1.0;
1669 real vrescale_energy(t_grpopts *opts, double therm_integral[])
1675 for (i = 0; i < opts->ngtc; i++)
1677 ener += therm_integral[i];
1683 void rescale_velocities(gmx_ekindata_t *ekind, t_mdatoms *mdatoms,
1684 int start, int end, rvec v[])
1687 t_grp_tcstat *tcstat;
1688 unsigned short *cACC, *cTC;
1693 tcstat = ekind->tcstat;
1698 gstat = ekind->grpstat;
1699 cACC = mdatoms->cACC;
1703 for (n = start; n < end; n++)
1713 /* Only scale the velocity component relative to the COM velocity */
1714 rvec_sub(v[n], gstat[ga].u, vrel);
1715 lg = tcstat[gt].lambda;
1716 for (d = 0; d < DIM; d++)
1718 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1725 for (n = start; n < end; n++)
1731 lg = tcstat[gt].lambda;
1732 for (d = 0; d < DIM; d++)
1741 /* set target temperatures if we are annealing */
1742 void update_annealing_target_temp(t_grpopts *opts, real t)
1744 int i, j, n, npoints;
1745 real pert, thist = 0, x;
1747 for (i = 0; i < opts->ngtc; i++)
1749 npoints = opts->anneal_npoints[i];
1750 switch (opts->annealing[i])
1755 /* calculate time modulo the period */
1756 pert = opts->anneal_time[i][npoints-1];
1758 thist = t - n*pert; /* modulo time */
1759 /* Make sure rounding didn't get us outside the interval */
1760 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1769 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)", i, opts->ngtc, npoints);
1771 /* We are doing annealing for this group if we got here,
1772 * and we have the (relative) time as thist.
1773 * calculate target temp */
1775 while ((j < npoints-1) && (thist > (opts->anneal_time[i][j+1])))
1781 /* Found our position between points j and j+1.
1782 * Interpolate: x is the amount from j+1, (1-x) from point j
1783 * First treat possible jumps in temperature as a special case.
1785 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1787 opts->ref_t[i] = opts->anneal_temp[i][j+1];
1791 x = ((thist-opts->anneal_time[i][j])/
1792 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1793 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1798 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];