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.
41 #include "gromacs/legacyheaders/typedefs.h"
42 #include "gromacs/legacyheaders/types/commrec.h"
43 #include "gromacs/utility/smalloc.h"
44 #include "gromacs/legacyheaders/update.h"
45 #include "gromacs/math/vec.h"
46 #include "gromacs/legacyheaders/macros.h"
47 #include "gromacs/math/units.h"
48 #include "gromacs/legacyheaders/names.h"
49 #include "gromacs/utility/fatalerror.h"
50 #include "gromacs/legacyheaders/txtdump.h"
51 #include "gromacs/legacyheaders/nrnb.h"
52 #include "gromacs/random/random.h"
53 #include "gromacs/legacyheaders/update.h"
54 #include "gromacs/legacyheaders/mdrun.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 void andersen_tcoupl(t_inputrec *ir, gmx_int64_t step,
726 const t_commrec *cr, const t_mdatoms *md, t_state *state, real rate, const gmx_bool *randomize, const real *boltzfac)
728 const int *gatindex = (DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
732 /* randomize the velocities of the selected particles */
734 for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */
736 int ng = gatindex ? gatindex[i] : i;
741 gc = md->cTC[i]; /* assign the atom to a temperature group if there are more than one */
745 if (ir->etc == etcANDERSEN)
753 gmx_rng_cycle_2uniform(step*2, ng, ir->andersen_seed, RND_SEED_ANDERSEN, uniform);
754 bRandomize = (uniform[0] < rate);
761 scal = sqrt(boltzfac[gc]*md->invmass[i]);
762 gmx_rng_cycle_3gaussian_table(step*2+1, ng, ir->andersen_seed, RND_SEED_ANDERSEN, gauss);
763 for (d = 0; d < DIM; d++)
765 state->v[i][d] = scal*gauss[d];
773 void nosehoover_tcoupl(t_grpopts *opts, gmx_ekindata_t *ekind, real dt,
774 double xi[], double vxi[], t_extmass *MassQ)
779 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
781 for (i = 0; (i < opts->ngtc); i++)
783 reft = max(0.0, opts->ref_t[i]);
785 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
786 xi[i] += dt*(oldvxi + vxi[i])*0.5;
790 t_state *init_bufstate(const t_state *template_state)
793 int nc = template_state->nhchainlength;
795 snew(state->nosehoover_xi, nc*template_state->ngtc);
796 snew(state->nosehoover_vxi, nc*template_state->ngtc);
797 snew(state->therm_integral, template_state->ngtc);
798 snew(state->nhpres_xi, nc*template_state->nnhpres);
799 snew(state->nhpres_vxi, nc*template_state->nnhpres);
804 void destroy_bufstate(t_state *state)
808 sfree(state->nosehoover_xi);
809 sfree(state->nosehoover_vxi);
810 sfree(state->therm_integral);
811 sfree(state->nhpres_xi);
812 sfree(state->nhpres_vxi);
816 void trotter_update(t_inputrec *ir, gmx_int64_t step, gmx_ekindata_t *ekind,
817 gmx_enerdata_t *enerd, t_state *state,
818 tensor vir, t_mdatoms *md,
819 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
822 int n, i, j, d, ntgrp, ngtc, gc = 0;
823 t_grp_tcstat *tcstat;
825 gmx_int64_t step_eff;
826 real ecorr, pcorr, dvdlcorr;
827 real bmass, qmass, reft, kT, dt, nd;
828 tensor dumpres, dumvir;
829 double *scalefac, dtc;
831 rvec sumv = {0, 0, 0}, consk;
834 if (trotter_seqno <= ettTSEQ2)
836 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
837 is actually the last half step from the previous step. Thus the first half step
838 actually corresponds to the n-1 step*/
846 bCouple = (ir->nsttcouple == 1 ||
847 do_per_step(step_eff+ir->nsttcouple, ir->nsttcouple));
849 trotter_seq = trotter_seqlist[trotter_seqno];
851 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
855 dtc = ir->nsttcouple*ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
856 opts = &(ir->opts); /* just for ease of referencing */
859 snew(scalefac, opts->ngtc);
860 for (i = 0; i < ngtc; i++)
864 /* execute the series of trotter updates specified in the trotterpart array */
866 for (i = 0; i < NTROTTERPARTS; i++)
868 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
869 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
878 switch (trotter_seq[i])
882 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir,
883 enerd->term[F_PDISPCORR], MassQ);
887 NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi,
888 state->nhpres_vxi, NULL, &(state->veta), MassQ, FALSE);
892 NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi,
893 state->nosehoover_vxi, scalefac, NULL, MassQ, (ir->eI == eiVV));
894 /* need to rescale the kinetic energies and velocities here. Could
895 scale the velocities later, but we need them scaled in order to
896 produce the correct outputs, so we'll scale them here. */
898 for (i = 0; i < ngtc; i++)
900 tcstat = &ekind->tcstat[i];
901 tcstat->vscale_nhc = scalefac[i];
902 tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
903 tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
905 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
906 /* but do we actually need the total? */
908 /* modify the velocities as well */
909 for (n = 0; n < md->homenr; n++)
911 if (md->cTC) /* does this conditional need to be here? is this always true?*/
915 for (d = 0; d < DIM; d++)
917 state->v[n][d] *= scalefac[gc];
922 for (d = 0; d < DIM; d++)
924 sumv[d] += (state->v[n][d])/md->invmass[n];
933 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
939 for (d = 0; d < DIM; d++)
941 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
943 fprintf(debug, "Conserved kappa: %15.8f %15.8f %15.8f\n", consk[0], consk[1], consk[2]);
951 extern void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit)
953 int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
954 t_grp_tcstat *tcstat;
956 real ecorr, pcorr, dvdlcorr;
957 real bmass, qmass, reft, kT, dt, ndj, nd;
958 tensor dumpres, dumvir;
960 opts = &(ir->opts); /* just for ease of referencing */
961 ngtc = ir->opts.ngtc;
962 nnhpres = state->nnhpres;
963 nh = state->nhchainlength;
969 snew(MassQ->Qinv, ngtc);
971 for (i = 0; (i < ngtc); i++)
973 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
975 MassQ->Qinv[i] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
979 MassQ->Qinv[i] = 0.0;
983 else if (EI_VV(ir->eI))
985 /* Set pressure variables */
989 if (state->vol0 == 0)
991 state->vol0 = det(state->box);
992 /* because we start by defining a fixed
993 compressibility, we need the volume at this
994 compressibility to solve the problem. */
998 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
999 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1000 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
1001 /* An alternate mass definition, from Tuckerman et al. */
1002 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1003 for (d = 0; d < DIM; d++)
1005 for (n = 0; n < DIM; n++)
1007 MassQ->Winvm[d][n] = PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
1008 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1009 before using MTTK for anisotropic states.*/
1012 /* Allocate space for thermostat variables */
1015 snew(MassQ->Qinv, ngtc*nh);
1018 /* now, set temperature variables */
1019 for (i = 0; i < ngtc; i++)
1021 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1023 reft = max(0.0, opts->ref_t[i]);
1026 for (j = 0; j < nh; j++)
1036 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
1042 for (j = 0; j < nh; j++)
1044 MassQ->Qinv[i*nh+j] = 0.0;
1051 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
1053 int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
1054 t_grp_tcstat *tcstat;
1056 real ecorr, pcorr, dvdlcorr;
1057 real bmass, qmass, reft, kT, dt, ndj, nd;
1058 tensor dumpres, dumvir;
1061 opts = &(ir->opts); /* just for ease of referencing */
1063 nnhpres = state->nnhpres;
1064 nh = state->nhchainlength;
1066 init_npt_masses(ir, state, MassQ, TRUE);
1068 /* first, initialize clear all the trotter calls */
1069 snew(trotter_seq, ettTSEQMAX);
1070 for (i = 0; i < ettTSEQMAX; i++)
1072 snew(trotter_seq[i], NTROTTERPARTS);
1073 for (j = 0; j < NTROTTERPARTS; j++)
1075 trotter_seq[i][j] = etrtNONE;
1077 trotter_seq[i][0] = etrtSKIPALL;
1082 /* no trotter calls, so we never use the values in the array.
1083 * We access them (so we need to define them, but ignore
1089 /* compute the kinetic energy by using the half step velocities or
1090 * the kinetic energies, depending on the order of the trotter calls */
1094 if (IR_NPT_TROTTER(ir))
1096 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1097 We start with the initial one. */
1098 /* first, a round that estimates veta. */
1099 trotter_seq[0][0] = etrtBAROV;
1101 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1103 /* The first half trotter update */
1104 trotter_seq[2][0] = etrtBAROV;
1105 trotter_seq[2][1] = etrtNHC;
1106 trotter_seq[2][2] = etrtBARONHC;
1108 /* The second half trotter update */
1109 trotter_seq[3][0] = etrtBARONHC;
1110 trotter_seq[3][1] = etrtNHC;
1111 trotter_seq[3][2] = etrtBAROV;
1113 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1116 else if (IR_NVT_TROTTER(ir))
1118 /* This is the easy version - there are only two calls, both the same.
1119 Otherwise, even easier -- no calls */
1120 trotter_seq[2][0] = etrtNHC;
1121 trotter_seq[3][0] = etrtNHC;
1123 else if (IR_NPH_TROTTER(ir))
1125 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1126 We start with the initial one. */
1127 /* first, a round that estimates veta. */
1128 trotter_seq[0][0] = etrtBAROV;
1130 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1132 /* The first half trotter update */
1133 trotter_seq[2][0] = etrtBAROV;
1134 trotter_seq[2][1] = etrtBARONHC;
1136 /* The second half trotter update */
1137 trotter_seq[3][0] = etrtBARONHC;
1138 trotter_seq[3][1] = etrtBAROV;
1140 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1143 else if (ir->eI == eiVVAK)
1145 if (IR_NPT_TROTTER(ir))
1147 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1148 We start with the initial one. */
1149 /* first, a round that estimates veta. */
1150 trotter_seq[0][0] = etrtBAROV;
1152 /* The first half trotter update, part 1 -- double update, because it commutes */
1153 trotter_seq[1][0] = etrtNHC;
1155 /* The first half trotter update, part 2 */
1156 trotter_seq[2][0] = etrtBAROV;
1157 trotter_seq[2][1] = etrtBARONHC;
1159 /* The second half trotter update, part 1 */
1160 trotter_seq[3][0] = etrtBARONHC;
1161 trotter_seq[3][1] = etrtBAROV;
1163 /* The second half trotter update */
1164 trotter_seq[4][0] = etrtNHC;
1166 else if (IR_NVT_TROTTER(ir))
1168 /* This is the easy version - there is only one call, both the same.
1169 Otherwise, even easier -- no calls */
1170 trotter_seq[1][0] = etrtNHC;
1171 trotter_seq[4][0] = etrtNHC;
1173 else if (IR_NPH_TROTTER(ir))
1175 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1176 We start with the initial one. */
1177 /* first, a round that estimates veta. */
1178 trotter_seq[0][0] = etrtBAROV;
1180 /* The first half trotter update, part 1 -- leave zero */
1181 trotter_seq[1][0] = etrtNHC;
1183 /* The first half trotter update, part 2 */
1184 trotter_seq[2][0] = etrtBAROV;
1185 trotter_seq[2][1] = etrtBARONHC;
1187 /* The second half trotter update, part 1 */
1188 trotter_seq[3][0] = etrtBARONHC;
1189 trotter_seq[3][1] = etrtBAROV;
1191 /* The second half trotter update -- blank for now */
1199 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
1202 snew(MassQ->QPinv, nnhpres*opts->nhchainlength);
1204 /* barostat temperature */
1205 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1207 reft = max(0.0, opts->ref_t[0]);
1209 for (i = 0; i < nnhpres; i++)
1211 for (j = 0; j < nh; j++)
1221 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1227 for (i = 0; i < nnhpres; i++)
1229 for (j = 0; j < nh; j++)
1231 MassQ->QPinv[i*nh+j] = 0.0;
1238 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1240 int i, j, nd, ndj, bmass, qmass, ngtcall;
1241 real ener_npt, reft, eta, kT, tau;
1244 real vol, dbaro, W, Q;
1245 int nh = state->nhchainlength;
1249 /* now we compute the contribution of the pressure to the conserved quantity*/
1251 if (ir->epc == epcMTTK)
1253 /* find the volume, and the kinetic energy of the volume */
1259 /* contribution from the pressure momenenta */
1260 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1262 /* contribution from the PV term */
1263 vol = det(state->box);
1264 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1267 case epctANISOTROPIC:
1271 case epctSURFACETENSION:
1274 case epctSEMIISOTROPIC:
1282 if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1284 /* add the energy from the barostat thermostat chain */
1285 for (i = 0; i < state->nnhpres; i++)
1288 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1289 ivxi = &state->nhpres_vxi[i*nh];
1290 ixi = &state->nhpres_xi[i*nh];
1291 iQinv = &(MassQ->QPinv[i*nh]);
1292 reft = max(ir->opts.ref_t[0], 0); /* using 'System' temperature */
1295 for (j = 0; j < nh; j++)
1299 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1300 /* contribution from the thermal variable of the NH chain */
1301 ener_npt += ixi[j]*kT;
1305 fprintf(debug, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i, j, ivxi[j], ixi[j]);
1313 for (i = 0; i < ir->opts.ngtc; i++)
1315 ixi = &state->nosehoover_xi[i*nh];
1316 ivxi = &state->nosehoover_vxi[i*nh];
1317 iQinv = &(MassQ->Qinv[i*nh]);
1319 nd = ir->opts.nrdf[i];
1320 reft = max(ir->opts.ref_t[i], 0);
1325 if (IR_NVT_TROTTER(ir))
1327 /* contribution from the thermal momenta of the NH chain */
1328 for (j = 0; j < nh; j++)
1332 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1333 /* contribution from the thermal variable of the NH chain */
1342 ener_npt += ndj*ixi[j]*kT;
1346 else /* Other non Trotter temperature NH control -- no chains yet. */
1348 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1349 ener_npt += nd*ixi[0]*kT;
1357 static real vrescale_gamdev(real ia,
1358 gmx_int64_t step, gmx_int64_t *count,
1359 gmx_int64_t seed1, gmx_int64_t seed2)
1360 /* Gamma distribution, adapted from numerical recipes */
1362 real am, e, s, v1, v2, x, y;
1373 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1375 v2 = 2.0*rnd[1] - 1.0;
1377 while (v1*v1 + v2*v2 > 1.0 ||
1378 v1*v1*GMX_REAL_MAX < 3.0*ia);
1379 /* The last check above ensures that both x (3.0 > 2.0 in s)
1380 * and the pre-factor for e do not go out of range.
1384 s = sqrt(2.0*am + 1.0);
1389 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1391 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1398 static real gaussian_count(gmx_int64_t step, gmx_int64_t *count,
1399 gmx_int64_t seed1, gmx_int64_t seed2)
1401 double rnd[2], x, y, r;
1405 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1406 x = 2.0*rnd[0] - 1.0;
1407 y = 2.0*rnd[1] - 1.0;
1410 while (r > 1.0 || r == 0.0);
1412 r = sqrt(-2.0*log(r)/r);
1417 static real vrescale_sumnoises(real nn,
1418 gmx_int64_t step, gmx_int64_t *count,
1419 gmx_int64_t seed1, gmx_int64_t seed2)
1422 * Returns the sum of nn independent gaussian noises squared
1423 * (i.e. equivalent to summing the square of the return values
1424 * of nn calls to gmx_rng_gaussian_real).
1426 const real ndeg_tol = 0.0001;
1429 if (nn < 2 + ndeg_tol)
1434 nn_int = (int)(nn + 0.5);
1436 if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1438 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);
1442 for (i = 0; i < nn_int; i++)
1444 gauss = gaussian_count(step, count, seed1, seed2);
1451 /* Use a gamma distribution for any real nn > 2 */
1452 r = 2.0*vrescale_gamdev(0.5*nn, step, count, seed1, seed2);
1458 static real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut,
1459 gmx_int64_t step, gmx_int64_t seed)
1462 * Generates a new value for the kinetic energy,
1463 * according to Bussi et al JCP (2007), Eq. (A7)
1464 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1465 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1466 * ndeg: number of degrees of freedom of the atoms to be thermalized
1467 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1469 /* rnd_count tracks the step-local state for the cycle random
1472 gmx_int64_t rnd_count = 0;
1473 real factor, rr, ekin_new;
1477 factor = exp(-1.0/taut);
1484 rr = gaussian_count(step, &rnd_count, seed, RND_SEED_VRESCALE);
1488 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, step, &rnd_count, seed, RND_SEED_VRESCALE) + rr*rr)/ndeg - kk) +
1489 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1494 void vrescale_tcoupl(t_inputrec *ir, gmx_int64_t step,
1495 gmx_ekindata_t *ekind, real dt,
1496 double therm_integral[])
1500 real Ek, Ek_ref1, Ek_ref, Ek_new;
1504 for (i = 0; (i < opts->ngtc); i++)
1508 Ek = trace(ekind->tcstat[i].ekinf);
1512 Ek = trace(ekind->tcstat[i].ekinh);
1515 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1517 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1518 Ek_ref = Ek_ref1*opts->nrdf[i];
1520 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i],
1524 /* Analytically Ek_new>=0, but we check for rounding errors */
1527 ekind->tcstat[i].lambda = 0.0;
1531 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1534 therm_integral[i] -= Ek_new - Ek;
1538 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1539 i, Ek_ref, Ek, Ek_new, ekind->tcstat[i].lambda);
1544 ekind->tcstat[i].lambda = 1.0;
1549 real vrescale_energy(t_grpopts *opts, double therm_integral[])
1555 for (i = 0; i < opts->ngtc; i++)
1557 ener += therm_integral[i];
1563 void rescale_velocities(gmx_ekindata_t *ekind, t_mdatoms *mdatoms,
1564 int start, int end, rvec v[])
1567 t_grp_tcstat *tcstat;
1568 unsigned short *cACC, *cTC;
1573 tcstat = ekind->tcstat;
1578 gstat = ekind->grpstat;
1579 cACC = mdatoms->cACC;
1583 for (n = start; n < end; n++)
1593 /* Only scale the velocity component relative to the COM velocity */
1594 rvec_sub(v[n], gstat[ga].u, vrel);
1595 lg = tcstat[gt].lambda;
1596 for (d = 0; d < DIM; d++)
1598 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1605 for (n = start; n < end; n++)
1611 lg = tcstat[gt].lambda;
1612 for (d = 0; d < DIM; d++)
1621 /* set target temperatures if we are annealing */
1622 void update_annealing_target_temp(t_grpopts *opts, real t)
1624 int i, j, n, npoints;
1625 real pert, thist = 0, x;
1627 for (i = 0; i < opts->ngtc; i++)
1629 npoints = opts->anneal_npoints[i];
1630 switch (opts->annealing[i])
1635 /* calculate time modulo the period */
1636 pert = opts->anneal_time[i][npoints-1];
1638 thist = t - n*pert; /* modulo time */
1639 /* Make sure rounding didn't get us outside the interval */
1640 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1649 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)", i, opts->ngtc, npoints);
1651 /* We are doing annealing for this group if we got here,
1652 * and we have the (relative) time as thist.
1653 * calculate target temp */
1655 while ((j < npoints-1) && (thist > (opts->anneal_time[i][j+1])))
1661 /* Found our position between points j and j+1.
1662 * Interpolate: x is the amount from j+1, (1-x) from point j
1663 * First treat possible jumps in temperature as a special case.
1665 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1667 opts->ref_t[i] = opts->anneal_temp[i][j+1];
1671 x = ((thist-opts->anneal_time[i][j])/
1672 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1673 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1678 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];