1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROwing Monsters And Cloning Shrimps
47 #include "gmx_fatal.h"
50 #include "gmx_random.h"
54 #define NTROTTERPARTS 3
56 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
58 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
59 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
61 #define MAX_SUZUKI_YOSHIDA_NUM 5
62 #define SUZUKI_YOSHIDA_NUM 5
64 static const double sy_const_1[] = { 1. };
65 static const double sy_const_3[] = { 0.828981543588751, -0.657963087177502, 0.828981543588751 };
66 static const double sy_const_5[] = { 0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065, 0.2967324292201065 };
68 static const double* sy_const[] = {
78 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
82 {0.828981543588751,-0.657963087177502,0.828981543588751},
84 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
87 /* these integration routines are only referenced inside this file */
88 static void NHC_trotter(t_grpopts *opts, int nvar, gmx_ekindata_t *ekind, real dtfull,
89 double xi[], double vxi[], double scalefac[], real *veta, t_extmass *MassQ, gmx_bool bEkinAveVel)
92 /* general routine for both barostat and thermostat nose hoover chains */
94 int i, j, mi, mj, jmax;
95 double Ekin, Efac, reft, kT, nd;
102 int mstepsi, mstepsj;
103 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
104 int nh = opts->nhchainlength;
107 mstepsi = mstepsj = ns;
109 /* if scalefac is NULL, we are doing the NHC of the barostat */
112 if (scalefac == NULL)
117 for (i = 0; i < nvar; i++)
120 /* make it easier to iterate by selecting
121 out the sub-array that corresponds to this T group */
127 iQinv = &(MassQ->QPinv[i*nh]);
128 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
129 reft = max(0.0, opts->ref_t[0]);
130 Ekin = sqr(*veta)/MassQ->Winv;
134 iQinv = &(MassQ->Qinv[i*nh]);
135 tcstat = &ekind->tcstat[i];
137 reft = max(0.0, opts->ref_t[i]);
140 Ekin = 2*trace(tcstat->ekinf)*tcstat->ekinscalef_nhc;
144 Ekin = 2*trace(tcstat->ekinh)*tcstat->ekinscaleh_nhc;
149 for (mi = 0; mi < mstepsi; mi++)
151 for (mj = 0; mj < mstepsj; mj++)
153 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
154 dt = sy_const[ns][mj] * dtfull / mstepsi;
156 /* compute the thermal forces */
157 GQ[0] = iQinv[0]*(Ekin - nd*kT);
159 for (j = 0; j < nh-1; j++)
163 /* we actually don't need to update here if we save the
164 state of the GQ, but it's easier to just recompute*/
165 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
173 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
174 for (j = nh-1; j > 0; j--)
176 Efac = exp(-0.125*dt*ivxi[j]);
177 ivxi[j-1] = Efac*(ivxi[j-1]*Efac + 0.25*dt*GQ[j-1]);
180 Efac = exp(-0.5*dt*ivxi[0]);
191 /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
192 able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
193 think about this a bit more . . . */
195 GQ[0] = iQinv[0]*(Ekin - nd*kT);
197 /* update thermostat positions */
198 for (j = 0; j < nh; j++)
200 ixi[j] += 0.5*dt*ivxi[j];
203 for (j = 0; j < nh-1; j++)
205 Efac = exp(-0.125*dt*ivxi[j+1]);
206 ivxi[j] = Efac*(ivxi[j]*Efac + 0.25*dt*GQ[j]);
209 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
216 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
223 static void boxv_trotter(t_inputrec *ir, real *veta, real dt, tensor box,
224 gmx_ekindata_t *ekind, tensor vir, real pcorr, t_extmass *MassQ)
229 int i, j, d, n, nwall;
231 tensor Winvm, ekinmod, localpres;
233 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
234 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
237 if (ir->epct == epctSEMIISOTROPIC)
246 /* eta is in pure units. veta is in units of ps^-1. GW is in
247 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
248 taken to use only RATIOS of eta in updating the volume. */
250 /* we take the partial pressure tensors, modify the
251 kinetic energy tensor, and recovert to pressure */
253 if (ir->opts.nrdf[0] == 0)
255 gmx_fatal(FARGS, "Barostat is coupled to a T-group with no degrees of freedom\n");
257 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
258 alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
259 alpha *= ekind->tcstat[0].ekinscalef_nhc;
260 msmul(ekind->ekin, alpha, ekinmod);
261 /* for now, we use Elr = 0, because if you want to get it right, you
262 really should be using PME. Maybe print a warning? */
264 pscal = calc_pres(ir->ePBC, nwall, box, ekinmod, vir, localpres)+pcorr;
267 GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
273 * This file implements temperature and pressure coupling algorithms:
274 * For now only the Weak coupling and the modified weak coupling.
276 * Furthermore computation of pressure and temperature is done here
280 real calc_pres(int ePBC, int nwall, matrix box, tensor ekin, tensor vir,
286 if (ePBC == epbcNONE || (ePBC == epbcXY && nwall != 2))
292 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
293 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
297 fac = PRESFAC*2.0/det(box);
298 for (n = 0; (n < DIM); n++)
300 for (m = 0; (m < DIM); m++)
302 pres[n][m] = (ekin[n][m] - vir[n][m])*fac;
308 pr_rvecs(debug, 0, "PC: pres", pres, DIM);
309 pr_rvecs(debug, 0, "PC: ekin", ekin, DIM);
310 pr_rvecs(debug, 0, "PC: vir ", vir, DIM);
311 pr_rvecs(debug, 0, "PC: box ", box, DIM);
314 return trace(pres)/DIM;
317 real calc_temp(real ekin, real nrdf)
321 return (2.0*ekin)/(nrdf*BOLTZ);
329 void parrinellorahman_pcoupl(FILE *fplog, gmx_large_int_t step,
330 t_inputrec *ir, real dt, tensor pres,
331 tensor box, tensor box_rel, tensor boxv,
332 tensor M, matrix mu, gmx_bool bFirstStep)
334 /* This doesn't do any coordinate updating. It just
335 * integrates the box vector equations from the calculated
336 * acceleration due to pressure difference. We also compute
337 * the tensor M which is used in update to couple the particle
338 * coordinates to the box vectors.
340 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
343 * M_nk = (h') * (h' * h + h' h) * h
345 * with the dots denoting time derivatives and h is the transformation from
346 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
347 * This also goes for the pressure and M tensors - they are transposed relative
348 * to ours. Our equation thus becomes:
351 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
353 * where b is the gromacs box matrix.
354 * Our box accelerations are given by
356 * b = vol/W inv(box') * (P-ref_P) (=h')
361 real vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
362 real atot, arel, change, maxchange, xy_pressure;
363 tensor invbox, pdiff, t1, t2;
367 m_inv_ur0(box, invbox);
371 /* Note that PRESFAC does not occur here.
372 * The pressure and compressibility always occur as a product,
373 * therefore the pressure unit drops out.
375 maxl = max(box[XX][XX], box[YY][YY]);
376 maxl = max(maxl, box[ZZ][ZZ]);
377 for (d = 0; d < DIM; d++)
379 for (n = 0; n < DIM; n++)
382 (4*M_PI*M_PI*ir->compress[d][n])/(3*ir->tau_p*ir->tau_p*maxl);
386 m_sub(pres, ir->ref_p, pdiff);
388 if (ir->epct == epctSURFACETENSION)
390 /* Unlike Berendsen coupling it might not be trivial to include a z
391 * pressure correction here? On the other hand we don't scale the
392 * box momentarily, but change accelerations, so it might not be crucial.
394 xy_pressure = 0.5*(pres[XX][XX]+pres[YY][YY]);
395 for (d = 0; d < ZZ; d++)
397 pdiff[d][d] = (xy_pressure-(pres[ZZ][ZZ]-ir->ref_p[d][d]/box[d][d]));
401 tmmul(invbox, pdiff, t1);
402 /* Move the off-diagonal elements of the 'force' to one side to ensure
403 * that we obey the box constraints.
405 for (d = 0; d < DIM; d++)
407 for (n = 0; n < d; n++)
409 t1[d][n] += t1[n][d];
416 case epctANISOTROPIC:
417 for (d = 0; d < DIM; d++)
419 for (n = 0; n <= d; n++)
421 t1[d][n] *= winv[d][n]*vol;
426 /* calculate total volume acceleration */
427 atot = box[XX][XX]*box[YY][YY]*t1[ZZ][ZZ]+
428 box[XX][XX]*t1[YY][YY]*box[ZZ][ZZ]+
429 t1[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
431 /* set all RELATIVE box accelerations equal, and maintain total V
433 for (d = 0; d < DIM; d++)
435 for (n = 0; n <= d; n++)
437 t1[d][n] = winv[0][0]*vol*arel*box[d][n];
441 case epctSEMIISOTROPIC:
442 case epctSURFACETENSION:
443 /* Note the correction to pdiff above for surftens. coupling */
445 /* calculate total XY volume acceleration */
446 atot = box[XX][XX]*t1[YY][YY]+t1[XX][XX]*box[YY][YY];
447 arel = atot/(2*box[XX][XX]*box[YY][YY]);
448 /* set RELATIVE XY box accelerations equal, and maintain total V
449 * change speed. Dont change the third box vector accelerations */
450 for (d = 0; d < ZZ; d++)
452 for (n = 0; n <= d; n++)
454 t1[d][n] = winv[d][n]*vol*arel*box[d][n];
457 for (n = 0; n < DIM; n++)
459 t1[ZZ][n] *= winv[d][n]*vol;
463 gmx_fatal(FARGS, "Parrinello-Rahman pressure coupling type %s "
464 "not supported yet\n", EPCOUPLTYPETYPE(ir->epct));
469 for (d = 0; d < DIM; d++)
471 for (n = 0; n <= d; n++)
473 boxv[d][n] += dt*t1[d][n];
475 /* We do NOT update the box vectors themselves here, since
476 * we need them for shifting later. It is instead done last
477 * in the update() routine.
480 /* Calculate the change relative to diagonal elements-
481 since it's perfectly ok for the off-diagonal ones to
482 be zero it doesn't make sense to check the change relative
486 change = fabs(dt*boxv[d][n]/box[d][d]);
488 if (change > maxchange)
495 if (maxchange > 0.01 && fplog)
499 "\nStep %s Warning: Pressure scaling more than 1%%. "
500 "This may mean your system\n is not yet equilibrated. "
501 "Use of Parrinello-Rahman pressure coupling during\n"
502 "equilibration can lead to simulation instability, "
503 "and is discouraged.\n",
504 gmx_step_str(step, buf));
508 preserve_box_shape(ir, box_rel, boxv);
510 mtmul(boxv, box, t1); /* t1=boxv * b' */
511 mmul(invbox, t1, t2);
512 mtmul(t2, invbox, M);
514 /* Determine the scaling matrix mu for the coordinates */
515 for (d = 0; d < DIM; d++)
517 for (n = 0; n <= d; n++)
519 t1[d][n] = box[d][n] + dt*boxv[d][n];
522 preserve_box_shape(ir, box_rel, t1);
523 /* t1 is the box at t+dt, determine mu as the relative change */
524 mmul_ur0(invbox, t1, mu);
527 void berendsen_pcoupl(FILE *fplog, gmx_large_int_t step,
528 t_inputrec *ir, real dt, tensor pres, matrix box,
532 real scalar_pressure, xy_pressure, p_corr_z;
533 char *ptr, buf[STRLEN];
536 * Calculate the scaling matrix mu
540 for (d = 0; d < DIM; d++)
542 scalar_pressure += pres[d][d]/DIM;
545 xy_pressure += pres[d][d]/(DIM-1);
548 /* Pressure is now in bar, everywhere. */
549 #define factor(d, m) (ir->compress[d][m]*dt/ir->tau_p)
551 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
552 * necessary for triclinic scaling
558 for (d = 0; d < DIM; d++)
560 mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
563 case epctSEMIISOTROPIC:
564 for (d = 0; d < ZZ; d++)
566 mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
569 1.0 - factor(ZZ, ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
571 case epctANISOTROPIC:
572 for (d = 0; d < DIM; d++)
574 for (n = 0; n < DIM; n++)
576 mu[d][n] = (d == n ? 1.0 : 0.0)
577 -factor(d, n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
581 case epctSURFACETENSION:
582 /* ir->ref_p[0/1] is the reference surface-tension times *
583 * the number of surfaces */
584 if (ir->compress[ZZ][ZZ])
586 p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
590 /* when the compressibity is zero, set the pressure correction *
591 * in the z-direction to zero to get the correct surface tension */
594 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
595 for (d = 0; d < DIM-1; d++)
597 mu[d][d] = 1.0 + factor(d, d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
598 - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
602 gmx_fatal(FARGS, "Berendsen pressure coupling type %s not supported yet\n",
603 EPCOUPLTYPETYPE(ir->epct));
606 /* To fullfill the orientation restrictions on triclinic boxes
607 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
608 * the other elements of mu to first order.
610 mu[YY][XX] += mu[XX][YY];
611 mu[ZZ][XX] += mu[XX][ZZ];
612 mu[ZZ][YY] += mu[YY][ZZ];
619 pr_rvecs(debug, 0, "PC: pres ", pres, 3);
620 pr_rvecs(debug, 0, "PC: mu ", mu, 3);
623 if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 ||
624 mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01 ||
625 mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
628 sprintf(buf, "\nStep %s Warning: pressure scaling more than 1%%, "
630 gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
633 fprintf(fplog, "%s", buf);
635 fprintf(stderr, "%s", buf);
639 void berendsen_pscale(t_inputrec *ir, matrix mu,
640 matrix box, matrix box_rel,
641 int start, int nr_atoms,
642 rvec x[], unsigned short cFREEZE[],
645 ivec *nFreeze = ir->opts.nFreeze;
648 /* Scale the positions */
649 for (n = start; n < start+nr_atoms; n++)
658 x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
662 x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
666 x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
669 /* compute final boxlengths */
670 for (d = 0; d < DIM; d++)
672 box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
673 box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
674 box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
677 preserve_box_shape(ir, box_rel, box);
679 /* (un)shifting should NOT be done after this,
680 * since the box vectors might have changed
682 inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
685 void berendsen_tcoupl(t_inputrec *ir, gmx_ekindata_t *ekind, real dt)
689 real T, reft = 0, lll;
693 for (i = 0; (i < opts->ngtc); i++)
697 T = ekind->tcstat[i].T;
701 T = ekind->tcstat[i].Th;
704 if ((opts->tau_t[i] > 0) && (T > 0.0))
706 reft = max(0.0, opts->ref_t[i]);
707 lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
708 ekind->tcstat[i].lambda = max(min(lll, 1.25), 0.8);
712 ekind->tcstat[i].lambda = 1.0;
717 fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n",
718 i, T, ekind->tcstat[i].lambda);
723 static int poisson_variate(real lambda, gmx_rng_t rng)
735 p *= gmx_rng_uniform_real(rng);
742 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)
745 int i, j, k, d, len, n, ngtc, gc = 0;
746 int nshake, nsettle, nrandom, nrand_group;
747 real boltz, scal, reft, prand;
750 /* convenience variables */
754 /* idef is only passed in if it's chance-per-particle andersen, so
755 it essentially serves as a boolean to determine which type of
756 andersen is being used */
760 /* randomly atoms to randomize. However, all constraint
761 groups have to have either all of the atoms or none of the
765 1. Select whether or not to randomize each atom to get the correct probability.
766 2. Cycle through the constraint groups.
767 2a. for each constraint group, determine the fraction f of that constraint group that are
768 chosen to be randomized.
769 2b. all atoms in the constraint group are randomized with probability f.
773 if ((rate < 0.05) && (md->homenr > 50))
775 /* if the rate is relatively high, use a standard method, if low rate,
777 /* poisson distributions approxmation, more efficient for
778 * low rates, fewer random numbers required */
779 nrandom = poisson_variate(md->homenr*rate, rng); /* how many do we randomize? Use Poisson. */
780 /* now we know how many, choose them randomly. No worries about repeats, at this rate, it's negligible.
781 worst thing that happens, it lowers the true rate an negligible amount */
782 for (i = 0; i < nrandom; i++)
784 randatom[(int)(gmx_rng_uniform_real(rng)*md->homenr)] = TRUE;
789 for (i = 0; i < md->homenr; i++)
791 if (gmx_rng_uniform_real(rng) < rate)
799 /* instead of looping over the constraint groups, if we had a
800 list of which atoms were in which constraint groups, we
801 could then loop over only the groups that are randomized
802 now. But that is not available now. Create later after
803 determining whether there actually is any slowing. */
805 /* first, loop through the settles to make sure all groups either entirely randomized, or not randomized. */
807 nsettle = idef->il[F_SETTLE].nr/2;
808 for (i = 0; i < nsettle; i++)
810 iatoms = idef->il[F_SETTLE].iatoms;
812 for (k = 0; k < 3; k++) /* settles are always 3 atoms, hardcoded */
814 if (randatom[iatoms[2*i+1]+k])
816 nrand_group++; /* count the number of atoms to be shaken in the settles group */
817 randatom[iatoms[2*i+1]+k] = FALSE;
823 prand = (nrand_group)/3.0; /* use this fraction to compute the probability the
824 whole group is randomized */
825 if (gmx_rng_uniform_real(rng) < prand)
827 for (k = 0; k < 3; k++)
829 randatom[iatoms[2*i+1]+k] = TRUE; /* mark them all to be randomized */
836 /* now loop through the shake groups */
838 for (i = 0; i < nshake; i++)
840 iatoms = &(idef->il[F_CONSTR].iatoms[sblock[i]]);
841 len = sblock[i+1]-sblock[i];
843 for (k = 0; k < len; k++)
846 { /* only 2/3 of the sblock items are atoms, the others are labels */
847 if (randatom[iatoms[k]])
850 randatom[iatoms[k]] = FALSE; /* need to mark it false here in case the atom is in more than
851 one group in the shake block */
858 prand = (nrand_group)/(1.0*(2*len/3));
859 if (gmx_rng_uniform_real(rng) < prand)
861 for (k = 0; k < len; k++)
864 { /* only 2/3 of the sblock items are atoms, the others are labels */
865 randatom[iatoms[k]] = TRUE; /* randomize all of them */
875 for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */
879 randatom_list[n] = i;
883 nrandom = n; /* there are some values of nrandom for which
884 this algorithm won't work; for example all
885 water molecules and nrandom =/= 3. Better to
886 recount and use this number (which we
887 calculate anyway: it will not affect
888 the average number of atoms accepted.
894 /* if it's andersen-massive, then randomize all the atoms */
895 nrandom = md->homenr;
896 for (i = 0; i < nrandom; i++)
898 randatom_list[i] = i;
902 /* randomize the velocities of the selected particles */
904 for (i = 0; i < nrandom; i++) /* now loop over the list of atoms */
906 n = randatom_list[i];
909 gc = md->cTC[n]; /* assign the atom to a temperature group if there are more than one */
913 scal = sqrt(boltzfac[gc]*md->invmass[n]);
914 for (d = 0; d < DIM; d++)
916 state->v[n][d] = scal*gmx_rng_gaussian_table(rng);
919 randatom[n] = FALSE; /* unmark this atom for randomization */
924 void nosehoover_tcoupl(t_grpopts *opts, gmx_ekindata_t *ekind, real dt,
925 double xi[], double vxi[], t_extmass *MassQ)
930 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
932 for (i = 0; (i < opts->ngtc); i++)
934 reft = max(0.0, opts->ref_t[i]);
936 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
937 xi[i] += dt*(oldvxi + vxi[i])*0.5;
941 t_state *init_bufstate(const t_state *template_state)
944 int nc = template_state->nhchainlength;
946 snew(state->nosehoover_xi, nc*template_state->ngtc);
947 snew(state->nosehoover_vxi, nc*template_state->ngtc);
948 snew(state->therm_integral, template_state->ngtc);
949 snew(state->nhpres_xi, nc*template_state->nnhpres);
950 snew(state->nhpres_vxi, nc*template_state->nnhpres);
955 void destroy_bufstate(t_state *state)
959 sfree(state->nosehoover_xi);
960 sfree(state->nosehoover_vxi);
961 sfree(state->therm_integral);
962 sfree(state->nhpres_xi);
963 sfree(state->nhpres_vxi);
967 void trotter_update(t_inputrec *ir, gmx_large_int_t step, gmx_ekindata_t *ekind,
968 gmx_enerdata_t *enerd, t_state *state,
969 tensor vir, t_mdatoms *md,
970 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
973 int n, i, j, d, ntgrp, ngtc, gc = 0;
974 t_grp_tcstat *tcstat;
976 gmx_large_int_t step_eff;
977 real ecorr, pcorr, dvdlcorr;
978 real bmass, qmass, reft, kT, dt, nd;
979 tensor dumpres, dumvir;
980 double *scalefac, dtc;
982 rvec sumv = {0, 0, 0}, consk;
985 if (trotter_seqno <= ettTSEQ2)
987 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
988 is actually the last half step from the previous step. Thus the first half step
989 actually corresponds to the n-1 step*/
997 bCouple = (ir->nsttcouple == 1 ||
998 do_per_step(step_eff+ir->nsttcouple, ir->nsttcouple));
1000 trotter_seq = trotter_seqlist[trotter_seqno];
1002 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
1006 dtc = ir->nsttcouple*ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
1007 opts = &(ir->opts); /* just for ease of referencing */
1010 snew(scalefac, opts->ngtc);
1011 for (i = 0; i < ngtc; i++)
1015 /* execute the series of trotter updates specified in the trotterpart array */
1017 for (i = 0; i < NTROTTERPARTS; i++)
1019 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
1020 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
1029 switch (trotter_seq[i])
1033 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir,
1034 enerd->term[F_PDISPCORR], MassQ);
1038 NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi,
1039 state->nhpres_vxi, NULL, &(state->veta), MassQ, FALSE);
1043 NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi,
1044 state->nosehoover_vxi, scalefac, NULL, MassQ, (ir->eI == eiVV));
1045 /* need to rescale the kinetic energies and velocities here. Could
1046 scale the velocities later, but we need them scaled in order to
1047 produce the correct outputs, so we'll scale them here. */
1049 for (i = 0; i < ngtc; i++)
1051 tcstat = &ekind->tcstat[i];
1052 tcstat->vscale_nhc = scalefac[i];
1053 tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
1054 tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
1056 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
1057 /* but do we actually need the total? */
1059 /* modify the velocities as well */
1060 for (n = md->start; n < md->start+md->homenr; n++)
1062 if (md->cTC) /* does this conditional need to be here? is this always true?*/
1066 for (d = 0; d < DIM; d++)
1068 state->v[n][d] *= scalefac[gc];
1073 for (d = 0; d < DIM; d++)
1075 sumv[d] += (state->v[n][d])/md->invmass[n];
1084 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
1090 for (d = 0; d < DIM; d++)
1092 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
1094 fprintf(debug, "Conserved kappa: %15.8f %15.8f %15.8f\n", consk[0], consk[1], consk[2]);
1102 extern void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit)
1104 int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
1105 t_grp_tcstat *tcstat;
1107 real ecorr, pcorr, dvdlcorr;
1108 real bmass, qmass, reft, kT, dt, ndj, nd;
1109 tensor dumpres, dumvir;
1111 opts = &(ir->opts); /* just for ease of referencing */
1112 ngtc = ir->opts.ngtc;
1113 nnhpres = state->nnhpres;
1114 nh = state->nhchainlength;
1120 snew(MassQ->Qinv, ngtc);
1122 for (i = 0; (i < ngtc); i++)
1124 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1126 MassQ->Qinv[i] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
1130 MassQ->Qinv[i] = 0.0;
1134 else if (EI_VV(ir->eI))
1136 /* Set pressure variables */
1140 if (state->vol0 == 0)
1142 state->vol0 = det(state->box);
1143 /* because we start by defining a fixed
1144 compressibility, we need the volume at this
1145 compressibility to solve the problem. */
1149 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1150 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1151 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
1152 /* An alternate mass definition, from Tuckerman et al. */
1153 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1154 for (d = 0; d < DIM; d++)
1156 for (n = 0; n < DIM; n++)
1158 MassQ->Winvm[d][n] = PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
1159 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1160 before using MTTK for anisotropic states.*/
1163 /* Allocate space for thermostat variables */
1166 snew(MassQ->Qinv, ngtc*nh);
1169 /* now, set temperature variables */
1170 for (i = 0; i < ngtc; i++)
1172 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1174 reft = max(0.0, opts->ref_t[i]);
1177 for (j = 0; j < nh; j++)
1187 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
1193 for (j = 0; j < nh; j++)
1195 MassQ->Qinv[i*nh+j] = 0.0;
1202 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
1204 int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
1205 t_grp_tcstat *tcstat;
1207 real ecorr, pcorr, dvdlcorr;
1208 real bmass, qmass, reft, kT, dt, ndj, nd;
1209 tensor dumpres, dumvir;
1212 opts = &(ir->opts); /* just for ease of referencing */
1214 nnhpres = state->nnhpres;
1215 nh = state->nhchainlength;
1217 init_npt_masses(ir, state, MassQ, TRUE);
1219 /* first, initialize clear all the trotter calls */
1220 snew(trotter_seq, ettTSEQMAX);
1221 for (i = 0; i < ettTSEQMAX; i++)
1223 snew(trotter_seq[i], NTROTTERPARTS);
1224 for (j = 0; j < NTROTTERPARTS; j++)
1226 trotter_seq[i][j] = etrtNONE;
1228 trotter_seq[i][0] = etrtSKIPALL;
1233 /* no trotter calls, so we never use the values in the array.
1234 * We access them (so we need to define them, but ignore
1240 /* compute the kinetic energy by using the half step velocities or
1241 * the kinetic energies, depending on the order of the trotter calls */
1245 if (IR_NPT_TROTTER(ir))
1247 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1248 We start with the initial one. */
1249 /* first, a round that estimates veta. */
1250 trotter_seq[0][0] = etrtBAROV;
1252 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1254 /* The first half trotter update */
1255 trotter_seq[2][0] = etrtBAROV;
1256 trotter_seq[2][1] = etrtNHC;
1257 trotter_seq[2][2] = etrtBARONHC;
1259 /* The second half trotter update */
1260 trotter_seq[3][0] = etrtBARONHC;
1261 trotter_seq[3][1] = etrtNHC;
1262 trotter_seq[3][2] = etrtBAROV;
1264 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1267 else if (IR_NVT_TROTTER(ir))
1269 /* This is the easy version - there are only two calls, both the same.
1270 Otherwise, even easier -- no calls */
1271 trotter_seq[2][0] = etrtNHC;
1272 trotter_seq[3][0] = etrtNHC;
1274 else if (IR_NPH_TROTTER(ir))
1276 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1277 We start with the initial one. */
1278 /* first, a round that estimates veta. */
1279 trotter_seq[0][0] = etrtBAROV;
1281 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1283 /* The first half trotter update */
1284 trotter_seq[2][0] = etrtBAROV;
1285 trotter_seq[2][1] = etrtBARONHC;
1287 /* The second half trotter update */
1288 trotter_seq[3][0] = etrtBARONHC;
1289 trotter_seq[3][1] = etrtBAROV;
1291 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1294 else if (ir->eI == eiVVAK)
1296 if (IR_NPT_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 -- double update, because it commutes */
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 */
1315 trotter_seq[4][0] = etrtNHC;
1317 else if (IR_NVT_TROTTER(ir))
1319 /* This is the easy version - there is only one call, both the same.
1320 Otherwise, even easier -- no calls */
1321 trotter_seq[1][0] = etrtNHC;
1322 trotter_seq[4][0] = etrtNHC;
1324 else if (IR_NPH_TROTTER(ir))
1326 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1327 We start with the initial one. */
1328 /* first, a round that estimates veta. */
1329 trotter_seq[0][0] = etrtBAROV;
1331 /* The first half trotter update, part 1 -- leave zero */
1332 trotter_seq[1][0] = etrtNHC;
1334 /* The first half trotter update, part 2 */
1335 trotter_seq[2][0] = etrtBAROV;
1336 trotter_seq[2][1] = etrtBARONHC;
1338 /* The second half trotter update, part 1 */
1339 trotter_seq[3][0] = etrtBARONHC;
1340 trotter_seq[3][1] = etrtBAROV;
1342 /* The second half trotter update -- blank for now */
1350 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
1353 snew(MassQ->QPinv, nnhpres*opts->nhchainlength);
1355 /* barostat temperature */
1356 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1358 reft = max(0.0, opts->ref_t[0]);
1360 for (i = 0; i < nnhpres; i++)
1362 for (j = 0; j < nh; j++)
1372 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1378 for (i = 0; i < nnhpres; i++)
1380 for (j = 0; j < nh; j++)
1382 MassQ->QPinv[i*nh+j] = 0.0;
1389 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1391 int i, j, nd, ndj, bmass, qmass, ngtcall;
1392 real ener_npt, reft, eta, kT, tau;
1395 real vol, dbaro, W, Q;
1396 int nh = state->nhchainlength;
1400 /* now we compute the contribution of the pressure to the conserved quantity*/
1402 if (ir->epc == epcMTTK)
1404 /* find the volume, and the kinetic energy of the volume */
1410 /* contribution from the pressure momenenta */
1411 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1413 /* contribution from the PV term */
1414 vol = det(state->box);
1415 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1418 case epctANISOTROPIC:
1422 case epctSURFACETENSION:
1425 case epctSEMIISOTROPIC:
1433 if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1435 /* add the energy from the barostat thermostat chain */
1436 for (i = 0; i < state->nnhpres; i++)
1439 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1440 ivxi = &state->nhpres_vxi[i*nh];
1441 ixi = &state->nhpres_xi[i*nh];
1442 iQinv = &(MassQ->QPinv[i*nh]);
1443 reft = max(ir->opts.ref_t[0], 0); /* using 'System' temperature */
1446 for (j = 0; j < nh; j++)
1450 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1451 /* contribution from the thermal variable of the NH chain */
1452 ener_npt += ixi[j]*kT;
1456 fprintf(debug, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i, j, ivxi[j], ixi[j]);
1464 for (i = 0; i < ir->opts.ngtc; i++)
1466 ixi = &state->nosehoover_xi[i*nh];
1467 ivxi = &state->nosehoover_vxi[i*nh];
1468 iQinv = &(MassQ->Qinv[i*nh]);
1470 nd = ir->opts.nrdf[i];
1471 reft = max(ir->opts.ref_t[i], 0);
1476 if (IR_NVT_TROTTER(ir))
1478 /* contribution from the thermal momenta of the NH chain */
1479 for (j = 0; j < nh; j++)
1483 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1484 /* contribution from the thermal variable of the NH chain */
1493 ener_npt += ndj*ixi[j]*kT;
1497 else /* Other non Trotter temperature NH control -- no chains yet. */
1499 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1500 ener_npt += nd*ixi[0]*kT;
1508 static real vrescale_gamdev(int ia, gmx_rng_t rng)
1509 /* Gamma distribution, adapted from numerical recipes */
1512 real am, e, s, v1, v2, x, y;
1519 for (j = 1; j <= ia; j++)
1521 x *= gmx_rng_uniform_real(rng);
1535 v1 = gmx_rng_uniform_real(rng);
1536 v2 = 2.0*gmx_rng_uniform_real(rng)-1.0;
1538 while (v1*v1 + v2*v2 > 1.0 ||
1539 v1*v1*GMX_REAL_MAX < 3.0*ia);
1540 /* The last check above ensures that both x (3.0 > 2.0 in s)
1541 * and the pre-factor for e do not go out of range.
1545 s = sqrt(2.0*am + 1.0);
1549 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1551 while (gmx_rng_uniform_real(rng) > e);
1557 static real vrescale_sumnoises(int nn, gmx_rng_t rng)
1560 * Returns the sum of n independent gaussian noises squared
1561 * (i.e. equivalent to summing the square of the return values
1562 * of nn calls to gmx_rng_gaussian_real).xs
1572 rr = gmx_rng_gaussian_real(rng);
1575 else if (nn % 2 == 0)
1577 return 2.0*vrescale_gamdev(nn/2, rng);
1581 rr = gmx_rng_gaussian_real(rng);
1582 return 2.0*vrescale_gamdev((nn-1)/2, rng) + rr*rr;
1586 static real vrescale_resamplekin(real kk, real sigma, int ndeg, real taut,
1590 * Generates a new value for the kinetic energy,
1591 * according to Bussi et al JCP (2007), Eq. (A7)
1592 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1593 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1594 * ndeg: number of degrees of freedom of the atoms to be thermalized
1595 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1601 factor = exp(-1.0/taut);
1607 rr = gmx_rng_gaussian_real(rng);
1610 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, rng) + rr*rr)/ndeg - kk) +
1611 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1614 void vrescale_tcoupl(t_inputrec *ir, gmx_ekindata_t *ekind, real dt,
1615 double therm_integral[], gmx_rng_t rng)
1619 real Ek, Ek_ref1, Ek_ref, Ek_new;
1623 for (i = 0; (i < opts->ngtc); i++)
1627 Ek = trace(ekind->tcstat[i].ekinf);
1631 Ek = trace(ekind->tcstat[i].ekinh);
1634 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1636 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1637 Ek_ref = Ek_ref1*opts->nrdf[i];
1639 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i],
1640 opts->tau_t[i]/dt, rng);
1642 /* Analytically Ek_new>=0, but we check for rounding errors */
1645 ekind->tcstat[i].lambda = 0.0;
1649 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1652 therm_integral[i] -= Ek_new - Ek;
1656 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1657 i, Ek_ref, Ek, Ek_new, ekind->tcstat[i].lambda);
1662 ekind->tcstat[i].lambda = 1.0;
1667 real vrescale_energy(t_grpopts *opts, double therm_integral[])
1673 for (i = 0; i < opts->ngtc; i++)
1675 ener += therm_integral[i];
1681 void rescale_velocities(gmx_ekindata_t *ekind, t_mdatoms *mdatoms,
1682 int start, int end, rvec v[])
1685 t_grp_tcstat *tcstat;
1686 unsigned short *cACC, *cTC;
1691 tcstat = ekind->tcstat;
1696 gstat = ekind->grpstat;
1697 cACC = mdatoms->cACC;
1701 for (n = start; n < end; n++)
1711 /* Only scale the velocity component relative to the COM velocity */
1712 rvec_sub(v[n], gstat[ga].u, vrel);
1713 lg = tcstat[gt].lambda;
1714 for (d = 0; d < DIM; d++)
1716 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1723 for (n = start; n < end; n++)
1729 lg = tcstat[gt].lambda;
1730 for (d = 0; d < DIM; d++)
1739 /* set target temperatures if we are annealing */
1740 void update_annealing_target_temp(t_grpopts *opts, real t)
1742 int i, j, n, npoints;
1743 real pert, thist = 0, x;
1745 for (i = 0; i < opts->ngtc; i++)
1747 npoints = opts->anneal_npoints[i];
1748 switch (opts->annealing[i])
1753 /* calculate time modulo the period */
1754 pert = opts->anneal_time[i][npoints-1];
1756 thist = t - n*pert; /* modulo time */
1757 /* Make sure rounding didn't get us outside the interval */
1758 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1767 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)", i, opts->ngtc, npoints);
1769 /* We are doing annealing for this group if we got here,
1770 * and we have the (relative) time as thist.
1771 * calculate target temp */
1773 while ((j < npoints-1) && (thist > (opts->anneal_time[i][j+1])))
1779 /* Found our position between points j and j+1.
1780 * Interpolate: x is the amount from j+1, (1-x) from point j
1781 * First treat possible jumps in temperature as a special case.
1783 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1785 opts->ref_t[i] = opts->anneal_temp[i][j+1];
1789 x = ((thist-opts->anneal_time[i][j])/
1790 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1791 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1796 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];