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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
50 #include "gmx_fatal.h"
53 #include "gmx_random.h"
57 #define NTROTTERPARTS 3
59 /* these integration routines are only referenced inside this file */
60 static void NHC_trotter(t_grpopts *opts, int nvar, gmx_ekindata_t *ekind, real dtfull,
61 double xi[], double vxi[], double scalefac[], real *veta, t_extmass *MassQ, gmx_bool bEkinAveVel)
64 /* general routine for both barostat and thermostat nose hoover chains */
66 int i, j, mi, mj, jmax;
67 double Ekin, Efac, reft, kT, nd;
75 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
76 int nh = opts->nhchainlength;
79 mstepsi = mstepsj = ns;
81 /* if scalefac is NULL, we are doing the NHC of the barostat */
89 for (i = 0; i < nvar; i++)
92 /* make it easier to iterate by selecting
93 out the sub-array that corresponds to this T group */
99 iQinv = &(MassQ->QPinv[i*nh]);
100 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
101 reft = max(0.0, opts->ref_t[0]);
102 Ekin = sqr(*veta)/MassQ->Winv;
106 iQinv = &(MassQ->Qinv[i*nh]);
107 tcstat = &ekind->tcstat[i];
109 reft = max(0.0, opts->ref_t[i]);
112 Ekin = 2*trace(tcstat->ekinf)*tcstat->ekinscalef_nhc;
116 Ekin = 2*trace(tcstat->ekinh)*tcstat->ekinscaleh_nhc;
121 for (mi = 0; mi < mstepsi; mi++)
123 for (mj = 0; mj < mstepsj; mj++)
125 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
126 dt = sy_const[ns][mj] * dtfull / mstepsi;
128 /* compute the thermal forces */
129 GQ[0] = iQinv[0]*(Ekin - nd*kT);
131 for (j = 0; j < nh-1; j++)
135 /* we actually don't need to update here if we save the
136 state of the GQ, but it's easier to just recompute*/
137 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
145 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
146 for (j = nh-1; j > 0; j--)
148 Efac = exp(-0.125*dt*ivxi[j]);
149 ivxi[j-1] = Efac*(ivxi[j-1]*Efac + 0.25*dt*GQ[j-1]);
152 Efac = exp(-0.5*dt*ivxi[0]);
163 /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
164 able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
165 think about this a bit more . . . */
167 GQ[0] = iQinv[0]*(Ekin - nd*kT);
169 /* update thermostat positions */
170 for (j = 0; j < nh; j++)
172 ixi[j] += 0.5*dt*ivxi[j];
175 for (j = 0; j < nh-1; j++)
177 Efac = exp(-0.125*dt*ivxi[j+1]);
178 ivxi[j] = Efac*(ivxi[j]*Efac + 0.25*dt*GQ[j]);
181 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
188 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
195 static void boxv_trotter(t_inputrec *ir, real *veta, real dt, tensor box,
196 gmx_ekindata_t *ekind, tensor vir, real pcorr, t_extmass *MassQ)
201 int i, j, d, n, nwall;
203 tensor Winvm, ekinmod, localpres;
205 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
206 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
209 if (ir->epct == epctSEMIISOTROPIC)
218 /* eta is in pure units. veta is in units of ps^-1. GW is in
219 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
220 taken to use only RATIOS of eta in updating the volume. */
222 /* we take the partial pressure tensors, modify the
223 kinetic energy tensor, and recovert to pressure */
225 if (ir->opts.nrdf[0] == 0)
227 gmx_fatal(FARGS, "Barostat is coupled to a T-group with no degrees of freedom\n");
229 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
230 alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
231 alpha *= ekind->tcstat[0].ekinscalef_nhc;
232 msmul(ekind->ekin, alpha, ekinmod);
233 /* for now, we use Elr = 0, because if you want to get it right, you
234 really should be using PME. Maybe print a warning? */
236 pscal = calc_pres(ir->ePBC, nwall, box, ekinmod, vir, localpres)+pcorr;
239 GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
245 * This file implements temperature and pressure coupling algorithms:
246 * For now only the Weak coupling and the modified weak coupling.
248 * Furthermore computation of pressure and temperature is done here
252 real calc_pres(int ePBC, int nwall, matrix box, tensor ekin, tensor vir,
258 if (ePBC == epbcNONE || (ePBC == epbcXY && nwall != 2))
264 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
265 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
269 fac = PRESFAC*2.0/det(box);
270 for (n = 0; (n < DIM); n++)
272 for (m = 0; (m < DIM); m++)
274 pres[n][m] = (ekin[n][m] - vir[n][m])*fac;
280 pr_rvecs(debug, 0, "PC: pres", pres, DIM);
281 pr_rvecs(debug, 0, "PC: ekin", ekin, DIM);
282 pr_rvecs(debug, 0, "PC: vir ", vir, DIM);
283 pr_rvecs(debug, 0, "PC: box ", box, DIM);
286 return trace(pres)/DIM;
289 real calc_temp(real ekin, real nrdf)
293 return (2.0*ekin)/(nrdf*BOLTZ);
301 void parrinellorahman_pcoupl(FILE *fplog, gmx_large_int_t step,
302 t_inputrec *ir, real dt, tensor pres,
303 tensor box, tensor box_rel, tensor boxv,
304 tensor M, matrix mu, gmx_bool bFirstStep)
306 /* This doesn't do any coordinate updating. It just
307 * integrates the box vector equations from the calculated
308 * acceleration due to pressure difference. We also compute
309 * the tensor M which is used in update to couple the particle
310 * coordinates to the box vectors.
312 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
315 * M_nk = (h') * (h' * h + h' h) * h
317 * with the dots denoting time derivatives and h is the transformation from
318 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
319 * This also goes for the pressure and M tensors - they are transposed relative
320 * to ours. Our equation thus becomes:
323 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
325 * where b is the gromacs box matrix.
326 * Our box accelerations are given by
328 * b = vol/W inv(box') * (P-ref_P) (=h')
333 real vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
334 real atot, arel, change, maxchange, xy_pressure;
335 tensor invbox, pdiff, t1, t2;
339 m_inv_ur0(box, invbox);
343 /* Note that PRESFAC does not occur here.
344 * The pressure and compressibility always occur as a product,
345 * therefore the pressure unit drops out.
347 maxl = max(box[XX][XX], box[YY][YY]);
348 maxl = max(maxl, box[ZZ][ZZ]);
349 for (d = 0; d < DIM; d++)
351 for (n = 0; n < DIM; n++)
354 (4*M_PI*M_PI*ir->compress[d][n])/(3*ir->tau_p*ir->tau_p*maxl);
358 m_sub(pres, ir->ref_p, pdiff);
360 if (ir->epct == epctSURFACETENSION)
362 /* Unlike Berendsen coupling it might not be trivial to include a z
363 * pressure correction here? On the other hand we don't scale the
364 * box momentarily, but change accelerations, so it might not be crucial.
366 xy_pressure = 0.5*(pres[XX][XX]+pres[YY][YY]);
367 for (d = 0; d < ZZ; d++)
369 pdiff[d][d] = (xy_pressure-(pres[ZZ][ZZ]-ir->ref_p[d][d]/box[d][d]));
373 tmmul(invbox, pdiff, t1);
374 /* Move the off-diagonal elements of the 'force' to one side to ensure
375 * that we obey the box constraints.
377 for (d = 0; d < DIM; d++)
379 for (n = 0; n < d; n++)
381 t1[d][n] += t1[n][d];
388 case epctANISOTROPIC:
389 for (d = 0; d < DIM; d++)
391 for (n = 0; n <= d; n++)
393 t1[d][n] *= winv[d][n]*vol;
398 /* calculate total volume acceleration */
399 atot = box[XX][XX]*box[YY][YY]*t1[ZZ][ZZ]+
400 box[XX][XX]*t1[YY][YY]*box[ZZ][ZZ]+
401 t1[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
403 /* set all RELATIVE box accelerations equal, and maintain total V
405 for (d = 0; d < DIM; d++)
407 for (n = 0; n <= d; n++)
409 t1[d][n] = winv[0][0]*vol*arel*box[d][n];
413 case epctSEMIISOTROPIC:
414 case epctSURFACETENSION:
415 /* Note the correction to pdiff above for surftens. coupling */
417 /* calculate total XY volume acceleration */
418 atot = box[XX][XX]*t1[YY][YY]+t1[XX][XX]*box[YY][YY];
419 arel = atot/(2*box[XX][XX]*box[YY][YY]);
420 /* set RELATIVE XY box accelerations equal, and maintain total V
421 * change speed. Dont change the third box vector accelerations */
422 for (d = 0; d < ZZ; d++)
424 for (n = 0; n <= d; n++)
426 t1[d][n] = winv[d][n]*vol*arel*box[d][n];
429 for (n = 0; n < DIM; n++)
431 t1[ZZ][n] *= winv[d][n]*vol;
435 gmx_fatal(FARGS, "Parrinello-Rahman pressure coupling type %s "
436 "not supported yet\n", EPCOUPLTYPETYPE(ir->epct));
441 for (d = 0; d < DIM; d++)
443 for (n = 0; n <= d; n++)
445 boxv[d][n] += dt*t1[d][n];
447 /* We do NOT update the box vectors themselves here, since
448 * we need them for shifting later. It is instead done last
449 * in the update() routine.
452 /* Calculate the change relative to diagonal elements-
453 since it's perfectly ok for the off-diagonal ones to
454 be zero it doesn't make sense to check the change relative
458 change = fabs(dt*boxv[d][n]/box[d][d]);
460 if (change > maxchange)
467 if (maxchange > 0.01 && fplog)
471 "\nStep %s Warning: Pressure scaling more than 1%%. "
472 "This may mean your system\n is not yet equilibrated. "
473 "Use of Parrinello-Rahman pressure coupling during\n"
474 "equilibration can lead to simulation instability, "
475 "and is discouraged.\n",
476 gmx_step_str(step, buf));
480 preserve_box_shape(ir, box_rel, boxv);
482 mtmul(boxv, box, t1); /* t1=boxv * b' */
483 mmul(invbox, t1, t2);
484 mtmul(t2, invbox, M);
486 /* Determine the scaling matrix mu for the coordinates */
487 for (d = 0; d < DIM; d++)
489 for (n = 0; n <= d; n++)
491 t1[d][n] = box[d][n] + dt*boxv[d][n];
494 preserve_box_shape(ir, box_rel, t1);
495 /* t1 is the box at t+dt, determine mu as the relative change */
496 mmul_ur0(invbox, t1, mu);
499 void berendsen_pcoupl(FILE *fplog, gmx_large_int_t step,
500 t_inputrec *ir, real dt, tensor pres, matrix box,
504 real scalar_pressure, xy_pressure, p_corr_z;
505 char *ptr, buf[STRLEN];
508 * Calculate the scaling matrix mu
512 for (d = 0; d < DIM; d++)
514 scalar_pressure += pres[d][d]/DIM;
517 xy_pressure += pres[d][d]/(DIM-1);
520 /* Pressure is now in bar, everywhere. */
521 #define factor(d, m) (ir->compress[d][m]*dt/ir->tau_p)
523 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
524 * necessary for triclinic scaling
530 for (d = 0; d < DIM; d++)
532 mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
535 case epctSEMIISOTROPIC:
536 for (d = 0; d < ZZ; d++)
538 mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
541 1.0 - factor(ZZ, ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
543 case epctANISOTROPIC:
544 for (d = 0; d < DIM; d++)
546 for (n = 0; n < DIM; n++)
548 mu[d][n] = (d == n ? 1.0 : 0.0)
549 -factor(d, n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
553 case epctSURFACETENSION:
554 /* ir->ref_p[0/1] is the reference surface-tension times *
555 * the number of surfaces */
556 if (ir->compress[ZZ][ZZ])
558 p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
562 /* when the compressibity is zero, set the pressure correction *
563 * in the z-direction to zero to get the correct surface tension */
566 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
567 for (d = 0; d < DIM-1; d++)
569 mu[d][d] = 1.0 + factor(d, d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
570 - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
574 gmx_fatal(FARGS, "Berendsen pressure coupling type %s not supported yet\n",
575 EPCOUPLTYPETYPE(ir->epct));
578 /* To fullfill the orientation restrictions on triclinic boxes
579 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
580 * the other elements of mu to first order.
582 mu[YY][XX] += mu[XX][YY];
583 mu[ZZ][XX] += mu[XX][ZZ];
584 mu[ZZ][YY] += mu[YY][ZZ];
591 pr_rvecs(debug, 0, "PC: pres ", pres, 3);
592 pr_rvecs(debug, 0, "PC: mu ", mu, 3);
595 if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 ||
596 mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01 ||
597 mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
600 sprintf(buf, "\nStep %s Warning: pressure scaling more than 1%%, "
602 gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
605 fprintf(fplog, "%s", buf);
607 fprintf(stderr, "%s", buf);
611 void berendsen_pscale(t_inputrec *ir, matrix mu,
612 matrix box, matrix box_rel,
613 int start, int nr_atoms,
614 rvec x[], unsigned short cFREEZE[],
617 ivec *nFreeze = ir->opts.nFreeze;
620 /* Scale the positions */
621 for (n = start; n < start+nr_atoms; n++)
630 x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
634 x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
638 x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
641 /* compute final boxlengths */
642 for (d = 0; d < DIM; d++)
644 box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
645 box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
646 box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
649 preserve_box_shape(ir, box_rel, box);
651 /* (un)shifting should NOT be done after this,
652 * since the box vectors might have changed
654 inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
657 void berendsen_tcoupl(t_inputrec *ir, gmx_ekindata_t *ekind, real dt)
661 real T, reft = 0, lll;
665 for (i = 0; (i < opts->ngtc); i++)
669 T = ekind->tcstat[i].T;
673 T = ekind->tcstat[i].Th;
676 if ((opts->tau_t[i] > 0) && (T > 0.0))
678 reft = max(0.0, opts->ref_t[i]);
679 lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
680 ekind->tcstat[i].lambda = max(min(lll, 1.25), 0.8);
684 ekind->tcstat[i].lambda = 1.0;
689 fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n",
690 i, T, ekind->tcstat[i].lambda);
695 static int poisson_variate(real lambda, gmx_rng_t rng)
707 p *= gmx_rng_uniform_real(rng);
714 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)
717 int i, j, k, d, len, n, ngtc, gc = 0;
718 int nshake, nsettle, nrandom, nrand_group;
719 real boltz, scal, reft, prand;
722 /* convenience variables */
726 /* idef is only passed in if it's chance-per-particle andersen, so
727 it essentially serves as a boolean to determine which type of
728 andersen is being used */
732 /* randomly atoms to randomize. However, all constraint
733 groups have to have either all of the atoms or none of the
737 1. Select whether or not to randomize each atom to get the correct probability.
738 2. Cycle through the constraint groups.
739 2a. for each constraint group, determine the fraction f of that constraint group that are
740 chosen to be randomized.
741 2b. all atoms in the constraint group are randomized with probability f.
745 if ((rate < 0.05) && (md->homenr > 50))
747 /* if the rate is relatively high, use a standard method, if low rate,
749 /* poisson distributions approxmation, more efficient for
750 * low rates, fewer random numbers required */
751 nrandom = poisson_variate(md->homenr*rate, rng); /* how many do we randomize? Use Poisson. */
752 /* now we know how many, choose them randomly. No worries about repeats, at this rate, it's negligible.
753 worst thing that happens, it lowers the true rate an negligible amount */
754 for (i = 0; i < nrandom; i++)
756 randatom[(int)(gmx_rng_uniform_real(rng)*md->homenr)] = TRUE;
761 for (i = 0; i < md->homenr; i++)
763 if (gmx_rng_uniform_real(rng) < rate)
771 /* instead of looping over the constraint groups, if we had a
772 list of which atoms were in which constraint groups, we
773 could then loop over only the groups that are randomized
774 now. But that is not available now. Create later after
775 determining whether there actually is any slowing. */
777 /* first, loop through the settles to make sure all groups either entirely randomized, or not randomized. */
779 nsettle = idef->il[F_SETTLE].nr/2;
780 for (i = 0; i < nsettle; i++)
782 iatoms = idef->il[F_SETTLE].iatoms;
784 for (k = 0; k < 3; k++) /* settles are always 3 atoms, hardcoded */
786 if (randatom[iatoms[2*i+1]+k])
788 nrand_group++; /* count the number of atoms to be shaken in the settles group */
789 randatom[iatoms[2*i+1]+k] = FALSE;
795 prand = (nrand_group)/3.0; /* use this fraction to compute the probability the
796 whole group is randomized */
797 if (gmx_rng_uniform_real(rng) < prand)
799 for (k = 0; k < 3; k++)
801 randatom[iatoms[2*i+1]+k] = TRUE; /* mark them all to be randomized */
808 /* now loop through the shake groups */
810 for (i = 0; i < nshake; i++)
812 iatoms = &(idef->il[F_CONSTR].iatoms[sblock[i]]);
813 len = sblock[i+1]-sblock[i];
815 for (k = 0; k < len; k++)
818 { /* only 2/3 of the sblock items are atoms, the others are labels */
819 if (randatom[iatoms[k]])
822 randatom[iatoms[k]] = FALSE; /* need to mark it false here in case the atom is in more than
823 one group in the shake block */
830 prand = (nrand_group)/(1.0*(2*len/3));
831 if (gmx_rng_uniform_real(rng) < prand)
833 for (k = 0; k < len; k++)
836 { /* only 2/3 of the sblock items are atoms, the others are labels */
837 randatom[iatoms[k]] = TRUE; /* randomize all of them */
847 for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */
851 randatom_list[n] = i;
855 nrandom = n; /* there are some values of nrandom for which
856 this algorithm won't work; for example all
857 water molecules and nrandom =/= 3. Better to
858 recount and use this number (which we
859 calculate anyway: it will not affect
860 the average number of atoms accepted.
866 /* if it's andersen-massive, then randomize all the atoms */
867 nrandom = md->homenr;
868 for (i = 0; i < nrandom; i++)
870 randatom_list[i] = i;
874 /* randomize the velocities of the selected particles */
876 for (i = 0; i < nrandom; i++) /* now loop over the list of atoms */
878 n = randatom_list[i];
881 gc = md->cTC[n]; /* assign the atom to a temperature group if there are more than one */
885 scal = sqrt(boltzfac[gc]*md->invmass[n]);
886 for (d = 0; d < DIM; d++)
888 state->v[n][d] = scal*gmx_rng_gaussian_table(rng);
891 randatom[n] = FALSE; /* unmark this atom for randomization */
896 void nosehoover_tcoupl(t_grpopts *opts, gmx_ekindata_t *ekind, real dt,
897 double xi[], double vxi[], t_extmass *MassQ)
902 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
904 for (i = 0; (i < opts->ngtc); i++)
906 reft = max(0.0, opts->ref_t[i]);
908 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
909 xi[i] += dt*(oldvxi + vxi[i])*0.5;
913 t_state *init_bufstate(const t_state *template_state)
916 int nc = template_state->nhchainlength;
918 snew(state->nosehoover_xi, nc*template_state->ngtc);
919 snew(state->nosehoover_vxi, nc*template_state->ngtc);
920 snew(state->therm_integral, template_state->ngtc);
921 snew(state->nhpres_xi, nc*template_state->nnhpres);
922 snew(state->nhpres_vxi, nc*template_state->nnhpres);
927 void destroy_bufstate(t_state *state)
931 sfree(state->nosehoover_xi);
932 sfree(state->nosehoover_vxi);
933 sfree(state->therm_integral);
934 sfree(state->nhpres_xi);
935 sfree(state->nhpres_vxi);
939 void trotter_update(t_inputrec *ir, gmx_large_int_t step, gmx_ekindata_t *ekind,
940 gmx_enerdata_t *enerd, t_state *state,
941 tensor vir, t_mdatoms *md,
942 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
945 int n, i, j, d, ntgrp, ngtc, gc = 0;
946 t_grp_tcstat *tcstat;
948 gmx_large_int_t step_eff;
949 real ecorr, pcorr, dvdlcorr;
950 real bmass, qmass, reft, kT, dt, nd;
951 tensor dumpres, dumvir;
952 double *scalefac, dtc;
954 rvec sumv = {0, 0, 0}, consk;
957 if (trotter_seqno <= ettTSEQ2)
959 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
960 is actually the last half step from the previous step. Thus the first half step
961 actually corresponds to the n-1 step*/
969 bCouple = (ir->nsttcouple == 1 ||
970 do_per_step(step_eff+ir->nsttcouple, ir->nsttcouple));
972 trotter_seq = trotter_seqlist[trotter_seqno];
974 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
978 dtc = ir->nsttcouple*ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
979 opts = &(ir->opts); /* just for ease of referencing */
982 snew(scalefac, opts->ngtc);
983 for (i = 0; i < ngtc; i++)
987 /* execute the series of trotter updates specified in the trotterpart array */
989 for (i = 0; i < NTROTTERPARTS; i++)
991 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
992 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
1001 switch (trotter_seq[i])
1005 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir,
1006 enerd->term[F_PDISPCORR], MassQ);
1010 NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi,
1011 state->nhpres_vxi, NULL, &(state->veta), MassQ, FALSE);
1015 NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi,
1016 state->nosehoover_vxi, scalefac, NULL, MassQ, (ir->eI == eiVV));
1017 /* need to rescale the kinetic energies and velocities here. Could
1018 scale the velocities later, but we need them scaled in order to
1019 produce the correct outputs, so we'll scale them here. */
1021 for (i = 0; i < ngtc; i++)
1023 tcstat = &ekind->tcstat[i];
1024 tcstat->vscale_nhc = scalefac[i];
1025 tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
1026 tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
1028 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
1029 /* but do we actually need the total? */
1031 /* modify the velocities as well */
1032 for (n = md->start; n < md->start+md->homenr; n++)
1034 if (md->cTC) /* does this conditional need to be here? is this always true?*/
1038 for (d = 0; d < DIM; d++)
1040 state->v[n][d] *= scalefac[gc];
1045 for (d = 0; d < DIM; d++)
1047 sumv[d] += (state->v[n][d])/md->invmass[n];
1056 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
1062 for (d = 0; d < DIM; d++)
1064 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
1066 fprintf(debug, "Conserved kappa: %15.8f %15.8f %15.8f\n", consk[0], consk[1], consk[2]);
1074 extern void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit)
1076 int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
1077 t_grp_tcstat *tcstat;
1079 real ecorr, pcorr, dvdlcorr;
1080 real bmass, qmass, reft, kT, dt, ndj, nd;
1081 tensor dumpres, dumvir;
1083 opts = &(ir->opts); /* just for ease of referencing */
1084 ngtc = ir->opts.ngtc;
1085 nnhpres = state->nnhpres;
1086 nh = state->nhchainlength;
1092 snew(MassQ->Qinv, ngtc);
1094 for (i = 0; (i < ngtc); i++)
1096 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1098 MassQ->Qinv[i] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
1102 MassQ->Qinv[i] = 0.0;
1106 else if (EI_VV(ir->eI))
1108 /* Set pressure variables */
1112 if (state->vol0 == 0)
1114 state->vol0 = det(state->box);
1115 /* because we start by defining a fixed
1116 compressibility, we need the volume at this
1117 compressibility to solve the problem. */
1121 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1122 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1123 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
1124 /* An alternate mass definition, from Tuckerman et al. */
1125 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1126 for (d = 0; d < DIM; d++)
1128 for (n = 0; n < DIM; n++)
1130 MassQ->Winvm[d][n] = PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
1131 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1132 before using MTTK for anisotropic states.*/
1135 /* Allocate space for thermostat variables */
1138 snew(MassQ->Qinv, ngtc*nh);
1141 /* now, set temperature variables */
1142 for (i = 0; i < ngtc; i++)
1144 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1146 reft = max(0.0, opts->ref_t[i]);
1149 for (j = 0; j < nh; j++)
1159 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
1165 for (j = 0; j < nh; j++)
1167 MassQ->Qinv[i*nh+j] = 0.0;
1174 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
1176 int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
1177 t_grp_tcstat *tcstat;
1179 real ecorr, pcorr, dvdlcorr;
1180 real bmass, qmass, reft, kT, dt, ndj, nd;
1181 tensor dumpres, dumvir;
1184 opts = &(ir->opts); /* just for ease of referencing */
1186 nnhpres = state->nnhpres;
1187 nh = state->nhchainlength;
1189 init_npt_masses(ir, state, MassQ, TRUE);
1191 /* first, initialize clear all the trotter calls */
1192 snew(trotter_seq, ettTSEQMAX);
1193 for (i = 0; i < ettTSEQMAX; i++)
1195 snew(trotter_seq[i], NTROTTERPARTS);
1196 for (j = 0; j < NTROTTERPARTS; j++)
1198 trotter_seq[i][j] = etrtNONE;
1200 trotter_seq[i][0] = etrtSKIPALL;
1205 /* no trotter calls, so we never use the values in the array.
1206 * We access them (so we need to define them, but ignore
1212 /* compute the kinetic energy by using the half step velocities or
1213 * the kinetic energies, depending on the order of the trotter calls */
1217 if (IR_NPT_TROTTER(ir))
1219 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1220 We start with the initial one. */
1221 /* first, a round that estimates veta. */
1222 trotter_seq[0][0] = etrtBAROV;
1224 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1226 /* The first half trotter update */
1227 trotter_seq[2][0] = etrtBAROV;
1228 trotter_seq[2][1] = etrtNHC;
1229 trotter_seq[2][2] = etrtBARONHC;
1231 /* The second half trotter update */
1232 trotter_seq[3][0] = etrtBARONHC;
1233 trotter_seq[3][1] = etrtNHC;
1234 trotter_seq[3][2] = etrtBAROV;
1236 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1239 else if (IR_NVT_TROTTER(ir))
1241 /* This is the easy version - there are only two calls, both the same.
1242 Otherwise, even easier -- no calls */
1243 trotter_seq[2][0] = etrtNHC;
1244 trotter_seq[3][0] = etrtNHC;
1246 else if (IR_NPH_TROTTER(ir))
1248 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1249 We start with the initial one. */
1250 /* first, a round that estimates veta. */
1251 trotter_seq[0][0] = etrtBAROV;
1253 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1255 /* The first half trotter update */
1256 trotter_seq[2][0] = etrtBAROV;
1257 trotter_seq[2][1] = etrtBARONHC;
1259 /* The second half trotter update */
1260 trotter_seq[3][0] = etrtBARONHC;
1261 trotter_seq[3][1] = etrtBAROV;
1263 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1266 else if (ir->eI == eiVVAK)
1268 if (IR_NPT_TROTTER(ir))
1270 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1271 We start with the initial one. */
1272 /* first, a round that estimates veta. */
1273 trotter_seq[0][0] = etrtBAROV;
1275 /* The first half trotter update, part 1 -- double update, because it commutes */
1276 trotter_seq[1][0] = etrtNHC;
1278 /* The first half trotter update, part 2 */
1279 trotter_seq[2][0] = etrtBAROV;
1280 trotter_seq[2][1] = etrtBARONHC;
1282 /* The second half trotter update, part 1 */
1283 trotter_seq[3][0] = etrtBARONHC;
1284 trotter_seq[3][1] = etrtBAROV;
1286 /* The second half trotter update */
1287 trotter_seq[4][0] = etrtNHC;
1289 else if (IR_NVT_TROTTER(ir))
1291 /* This is the easy version - there is only one call, both the same.
1292 Otherwise, even easier -- no calls */
1293 trotter_seq[1][0] = etrtNHC;
1294 trotter_seq[4][0] = etrtNHC;
1296 else if (IR_NPH_TROTTER(ir))
1298 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1299 We start with the initial one. */
1300 /* first, a round that estimates veta. */
1301 trotter_seq[0][0] = etrtBAROV;
1303 /* The first half trotter update, part 1 -- leave zero */
1304 trotter_seq[1][0] = etrtNHC;
1306 /* The first half trotter update, part 2 */
1307 trotter_seq[2][0] = etrtBAROV;
1308 trotter_seq[2][1] = etrtBARONHC;
1310 /* The second half trotter update, part 1 */
1311 trotter_seq[3][0] = etrtBARONHC;
1312 trotter_seq[3][1] = etrtBAROV;
1314 /* The second half trotter update -- blank for now */
1322 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
1325 snew(MassQ->QPinv, nnhpres*opts->nhchainlength);
1327 /* barostat temperature */
1328 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1330 reft = max(0.0, opts->ref_t[0]);
1332 for (i = 0; i < nnhpres; i++)
1334 for (j = 0; j < nh; j++)
1344 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1350 for (i = 0; i < nnhpres; i++)
1352 for (j = 0; j < nh; j++)
1354 MassQ->QPinv[i*nh+j] = 0.0;
1361 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1363 int i, j, nd, ndj, bmass, qmass, ngtcall;
1364 real ener_npt, reft, eta, kT, tau;
1367 real vol, dbaro, W, Q;
1368 int nh = state->nhchainlength;
1372 /* now we compute the contribution of the pressure to the conserved quantity*/
1374 if (ir->epc == epcMTTK)
1376 /* find the volume, and the kinetic energy of the volume */
1382 /* contribution from the pressure momenenta */
1383 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1385 /* contribution from the PV term */
1386 vol = det(state->box);
1387 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1390 case epctANISOTROPIC:
1394 case epctSURFACETENSION:
1397 case epctSEMIISOTROPIC:
1405 if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1407 /* add the energy from the barostat thermostat chain */
1408 for (i = 0; i < state->nnhpres; i++)
1411 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1412 ivxi = &state->nhpres_vxi[i*nh];
1413 ixi = &state->nhpres_xi[i*nh];
1414 iQinv = &(MassQ->QPinv[i*nh]);
1415 reft = max(ir->opts.ref_t[0], 0); /* using 'System' temperature */
1418 for (j = 0; j < nh; j++)
1422 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1423 /* contribution from the thermal variable of the NH chain */
1424 ener_npt += ixi[j]*kT;
1428 fprintf(debug, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i, j, ivxi[j], ixi[j]);
1436 for (i = 0; i < ir->opts.ngtc; i++)
1438 ixi = &state->nosehoover_xi[i*nh];
1439 ivxi = &state->nosehoover_vxi[i*nh];
1440 iQinv = &(MassQ->Qinv[i*nh]);
1442 nd = ir->opts.nrdf[i];
1443 reft = max(ir->opts.ref_t[i], 0);
1448 if (IR_NVT_TROTTER(ir))
1450 /* contribution from the thermal momenta of the NH chain */
1451 for (j = 0; j < nh; j++)
1455 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1456 /* contribution from the thermal variable of the NH chain */
1465 ener_npt += ndj*ixi[j]*kT;
1469 else /* Other non Trotter temperature NH control -- no chains yet. */
1471 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1472 ener_npt += nd*ixi[0]*kT;
1480 static real vrescale_gamdev(real ia, gmx_rng_t rng)
1481 /* Gamma distribution, adapted from numerical recipes */
1483 real am, e, s, v1, v2, x, y;
1493 v1 = gmx_rng_uniform_real(rng);
1494 v2 = 2.0*gmx_rng_uniform_real(rng) - 1.0;
1496 while (v1*v1 + v2*v2 > 1.0 ||
1497 v1*v1*GMX_REAL_MAX < 3.0*ia);
1498 /* The last check above ensures that both x (3.0 > 2.0 in s)
1499 * and the pre-factor for e do not go out of range.
1503 s = sqrt(2.0*am + 1.0);
1508 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1510 while (gmx_rng_uniform_real(rng) > e);
1515 static real vrescale_sumnoises(real nn, gmx_rng_t rng)
1518 * Returns the sum of n independent gaussian noises squared
1519 * (i.e. equivalent to summing the square of the return values
1520 * of nn calls to gmx_rng_gaussian_real).xs
1522 const real ndeg_tol = 0.0001;
1525 if (nn < 2 + ndeg_tol)
1530 nn_int = (int)(nn + 0.5);
1532 if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1534 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);
1538 for (i = 0; i < nn_int; i++)
1540 gauss = gmx_rng_gaussian_real(rng);
1547 /* This call to the rng is only here to keep the rng state in sync
1548 * with pre-4.6.6 versions, so simulations results will be similar.
1549 * Remove this for 5.0.
1551 if (((int)nn) % 2 == 1)
1553 gmx_rng_gaussian_real(rng);
1556 /* Use a gamma distribution for any real nn > 2 */
1557 r = 2.0*vrescale_gamdev(0.5*nn, rng);
1563 static real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut,
1567 * Generates a new value for the kinetic energy,
1568 * according to Bussi et al JCP (2007), Eq. (A7)
1569 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1570 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1571 * ndeg: number of degrees of freedom of the atoms to be thermalized
1572 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1578 factor = exp(-1.0/taut);
1584 rr = gmx_rng_gaussian_real(rng);
1587 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, rng) + rr*rr)/ndeg - kk) +
1588 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1591 void vrescale_tcoupl(t_inputrec *ir, gmx_ekindata_t *ekind, real dt,
1592 double therm_integral[], gmx_rng_t rng)
1596 real Ek, Ek_ref1, Ek_ref, Ek_new;
1600 for (i = 0; (i < opts->ngtc); i++)
1604 Ek = trace(ekind->tcstat[i].ekinf);
1608 Ek = trace(ekind->tcstat[i].ekinh);
1611 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1613 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1614 Ek_ref = Ek_ref1*opts->nrdf[i];
1616 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i],
1617 opts->tau_t[i]/dt, rng);
1619 /* Analytically Ek_new>=0, but we check for rounding errors */
1622 ekind->tcstat[i].lambda = 0.0;
1626 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1629 therm_integral[i] -= Ek_new - Ek;
1633 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1634 i, Ek_ref, Ek, Ek_new, ekind->tcstat[i].lambda);
1639 ekind->tcstat[i].lambda = 1.0;
1644 real vrescale_energy(t_grpopts *opts, double therm_integral[])
1650 for (i = 0; i < opts->ngtc; i++)
1652 ener += therm_integral[i];
1658 void rescale_velocities(gmx_ekindata_t *ekind, t_mdatoms *mdatoms,
1659 int start, int end, rvec v[])
1662 t_grp_tcstat *tcstat;
1663 unsigned short *cACC, *cTC;
1668 tcstat = ekind->tcstat;
1673 gstat = ekind->grpstat;
1674 cACC = mdatoms->cACC;
1678 for (n = start; n < end; n++)
1688 /* Only scale the velocity component relative to the COM velocity */
1689 rvec_sub(v[n], gstat[ga].u, vrel);
1690 lg = tcstat[gt].lambda;
1691 for (d = 0; d < DIM; d++)
1693 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1700 for (n = start; n < end; n++)
1706 lg = tcstat[gt].lambda;
1707 for (d = 0; d < DIM; d++)
1716 /* set target temperatures if we are annealing */
1717 void update_annealing_target_temp(t_grpopts *opts, real t)
1719 int i, j, n, npoints;
1720 real pert, thist = 0, x;
1722 for (i = 0; i < opts->ngtc; i++)
1724 npoints = opts->anneal_npoints[i];
1725 switch (opts->annealing[i])
1730 /* calculate time modulo the period */
1731 pert = opts->anneal_time[i][npoints-1];
1733 thist = t - n*pert; /* modulo time */
1734 /* Make sure rounding didn't get us outside the interval */
1735 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1744 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)", i, opts->ngtc, npoints);
1746 /* We are doing annealing for this group if we got here,
1747 * and we have the (relative) time as thist.
1748 * calculate target temp */
1750 while ((j < npoints-1) && (thist > (opts->anneal_time[i][j+1])))
1756 /* Found our position between points j and j+1.
1757 * Interpolate: x is the amount from j+1, (1-x) from point j
1758 * First treat possible jumps in temperature as a special case.
1760 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1762 opts->ref_t[i] = opts->anneal_temp[i][j+1];
1766 x = ((thist-opts->anneal_time[i][j])/
1767 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1768 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1773 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];