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/macros.h"
42 #include "gromacs/legacyheaders/mdrun.h"
43 #include "gromacs/legacyheaders/names.h"
44 #include "gromacs/legacyheaders/nrnb.h"
45 #include "gromacs/legacyheaders/txtdump.h"
46 #include "gromacs/legacyheaders/typedefs.h"
47 #include "gromacs/legacyheaders/update.h"
48 #include "gromacs/legacyheaders/types/commrec.h"
49 #include "gromacs/math/units.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/random/random.h"
52 #include "gromacs/utility/fatalerror.h"
53 #include "gromacs/utility/smalloc.h"
55 #define NTROTTERPARTS 3
57 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
59 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
60 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
62 #define MAX_SUZUKI_YOSHIDA_NUM 5
63 #define SUZUKI_YOSHIDA_NUM 5
65 static const double sy_const_1[] = { 1. };
66 static const double sy_const_3[] = { 0.828981543588751, -0.657963087177502, 0.828981543588751 };
67 static const double sy_const_5[] = { 0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065, 0.2967324292201065 };
69 static const double* sy_const[] = {
79 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
83 {0.828981543588751,-0.657963087177502,0.828981543588751},
85 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
88 /* these integration routines are only referenced inside this file */
89 static void NHC_trotter(t_grpopts *opts, int nvar, gmx_ekindata_t *ekind, real dtfull,
90 double xi[], double vxi[], double scalefac[], real *veta, t_extmass *MassQ, gmx_bool bEkinAveVel)
93 /* general routine for both barostat and thermostat nose hoover chains */
95 int i, j, mi, mj, jmax;
96 double Ekin, Efac, reft, kT, nd;
103 int mstepsi, mstepsj;
104 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
105 int nh = opts->nhchainlength;
108 mstepsi = mstepsj = ns;
110 /* if scalefac is NULL, we are doing the NHC of the barostat */
113 if (scalefac == NULL)
118 for (i = 0; i < nvar; i++)
121 /* make it easier to iterate by selecting
122 out the sub-array that corresponds to this T group */
128 iQinv = &(MassQ->QPinv[i*nh]);
129 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
130 reft = max(0.0, opts->ref_t[0]);
131 Ekin = sqr(*veta)/MassQ->Winv;
135 iQinv = &(MassQ->Qinv[i*nh]);
136 tcstat = &ekind->tcstat[i];
138 reft = max(0.0, opts->ref_t[i]);
141 Ekin = 2*trace(tcstat->ekinf)*tcstat->ekinscalef_nhc;
145 Ekin = 2*trace(tcstat->ekinh)*tcstat->ekinscaleh_nhc;
150 for (mi = 0; mi < mstepsi; mi++)
152 for (mj = 0; mj < mstepsj; mj++)
154 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
155 dt = sy_const[ns][mj] * dtfull / mstepsi;
157 /* compute the thermal forces */
158 GQ[0] = iQinv[0]*(Ekin - nd*kT);
160 for (j = 0; j < nh-1; j++)
164 /* we actually don't need to update here if we save the
165 state of the GQ, but it's easier to just recompute*/
166 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
174 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
175 for (j = nh-1; j > 0; j--)
177 Efac = exp(-0.125*dt*ivxi[j]);
178 ivxi[j-1] = Efac*(ivxi[j-1]*Efac + 0.25*dt*GQ[j-1]);
181 Efac = exp(-0.5*dt*ivxi[0]);
192 /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
193 able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
194 think about this a bit more . . . */
196 GQ[0] = iQinv[0]*(Ekin - nd*kT);
198 /* update thermostat positions */
199 for (j = 0; j < nh; j++)
201 ixi[j] += 0.5*dt*ivxi[j];
204 for (j = 0; j < nh-1; j++)
206 Efac = exp(-0.125*dt*ivxi[j+1]);
207 ivxi[j] = Efac*(ivxi[j]*Efac + 0.25*dt*GQ[j]);
210 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
217 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
224 static void boxv_trotter(t_inputrec *ir, real *veta, real dt, tensor box,
225 gmx_ekindata_t *ekind, tensor vir, real pcorr, t_extmass *MassQ)
230 int i, j, d, n, nwall;
232 tensor Winvm, ekinmod, localpres;
234 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
235 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
238 if (ir->epct == epctSEMIISOTROPIC)
247 /* eta is in pure units. veta is in units of ps^-1. GW is in
248 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
249 taken to use only RATIOS of eta in updating the volume. */
251 /* we take the partial pressure tensors, modify the
252 kinetic energy tensor, and recovert to pressure */
254 if (ir->opts.nrdf[0] == 0)
256 gmx_fatal(FARGS, "Barostat is coupled to a T-group with no degrees of freedom\n");
258 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
259 alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
260 alpha *= ekind->tcstat[0].ekinscalef_nhc;
261 msmul(ekind->ekin, alpha, ekinmod);
262 /* for now, we use Elr = 0, because if you want to get it right, you
263 really should be using PME. Maybe print a warning? */
265 pscal = calc_pres(ir->ePBC, nwall, box, ekinmod, vir, localpres)+pcorr;
268 GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
274 * This file implements temperature and pressure coupling algorithms:
275 * For now only the Weak coupling and the modified weak coupling.
277 * Furthermore computation of pressure and temperature is done here
281 real calc_pres(int ePBC, int nwall, matrix box, tensor ekin, tensor vir,
287 if (ePBC == epbcNONE || (ePBC == epbcXY && nwall != 2))
293 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
294 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
298 fac = PRESFAC*2.0/det(box);
299 for (n = 0; (n < DIM); n++)
301 for (m = 0; (m < DIM); m++)
303 pres[n][m] = (ekin[n][m] - vir[n][m])*fac;
309 pr_rvecs(debug, 0, "PC: pres", pres, DIM);
310 pr_rvecs(debug, 0, "PC: ekin", ekin, DIM);
311 pr_rvecs(debug, 0, "PC: vir ", vir, DIM);
312 pr_rvecs(debug, 0, "PC: box ", box, DIM);
315 return trace(pres)/DIM;
318 real calc_temp(real ekin, real nrdf)
322 return (2.0*ekin)/(nrdf*BOLTZ);
330 void parrinellorahman_pcoupl(FILE *fplog, gmx_int64_t step,
331 t_inputrec *ir, real dt, tensor pres,
332 tensor box, tensor box_rel, tensor boxv,
333 tensor M, matrix mu, gmx_bool bFirstStep)
335 /* This doesn't do any coordinate updating. It just
336 * integrates the box vector equations from the calculated
337 * acceleration due to pressure difference. We also compute
338 * the tensor M which is used in update to couple the particle
339 * coordinates to the box vectors.
341 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
344 * M_nk = (h') * (h' * h + h' h) * h
346 * with the dots denoting time derivatives and h is the transformation from
347 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
348 * This also goes for the pressure and M tensors - they are transposed relative
349 * to ours. Our equation thus becomes:
352 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
354 * where b is the gromacs box matrix.
355 * Our box accelerations are given by
357 * b = vol/W inv(box') * (P-ref_P) (=h')
362 real vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
363 real atot, arel, change, maxchange, xy_pressure;
364 tensor invbox, pdiff, t1, t2;
368 m_inv_ur0(box, invbox);
372 /* Note that PRESFAC does not occur here.
373 * The pressure and compressibility always occur as a product,
374 * therefore the pressure unit drops out.
376 maxl = max(box[XX][XX], box[YY][YY]);
377 maxl = max(maxl, box[ZZ][ZZ]);
378 for (d = 0; d < DIM; d++)
380 for (n = 0; n < DIM; n++)
383 (4*M_PI*M_PI*ir->compress[d][n])/(3*ir->tau_p*ir->tau_p*maxl);
387 m_sub(pres, ir->ref_p, pdiff);
389 if (ir->epct == epctSURFACETENSION)
391 /* Unlike Berendsen coupling it might not be trivial to include a z
392 * pressure correction here? On the other hand we don't scale the
393 * box momentarily, but change accelerations, so it might not be crucial.
395 xy_pressure = 0.5*(pres[XX][XX]+pres[YY][YY]);
396 for (d = 0; d < ZZ; d++)
398 pdiff[d][d] = (xy_pressure-(pres[ZZ][ZZ]-ir->ref_p[d][d]/box[d][d]));
402 tmmul(invbox, pdiff, t1);
403 /* Move the off-diagonal elements of the 'force' to one side to ensure
404 * that we obey the box constraints.
406 for (d = 0; d < DIM; d++)
408 for (n = 0; n < d; n++)
410 t1[d][n] += t1[n][d];
417 case epctANISOTROPIC:
418 for (d = 0; d < DIM; d++)
420 for (n = 0; n <= d; n++)
422 t1[d][n] *= winv[d][n]*vol;
427 /* calculate total volume acceleration */
428 atot = box[XX][XX]*box[YY][YY]*t1[ZZ][ZZ]+
429 box[XX][XX]*t1[YY][YY]*box[ZZ][ZZ]+
430 t1[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
432 /* set all RELATIVE box accelerations equal, and maintain total V
434 for (d = 0; d < DIM; d++)
436 for (n = 0; n <= d; n++)
438 t1[d][n] = winv[0][0]*vol*arel*box[d][n];
442 case epctSEMIISOTROPIC:
443 case epctSURFACETENSION:
444 /* Note the correction to pdiff above for surftens. coupling */
446 /* calculate total XY volume acceleration */
447 atot = box[XX][XX]*t1[YY][YY]+t1[XX][XX]*box[YY][YY];
448 arel = atot/(2*box[XX][XX]*box[YY][YY]);
449 /* set RELATIVE XY box accelerations equal, and maintain total V
450 * change speed. Dont change the third box vector accelerations */
451 for (d = 0; d < ZZ; d++)
453 for (n = 0; n <= d; n++)
455 t1[d][n] = winv[d][n]*vol*arel*box[d][n];
458 for (n = 0; n < DIM; n++)
460 t1[ZZ][n] *= winv[d][n]*vol;
464 gmx_fatal(FARGS, "Parrinello-Rahman pressure coupling type %s "
465 "not supported yet\n", EPCOUPLTYPETYPE(ir->epct));
470 for (d = 0; d < DIM; d++)
472 for (n = 0; n <= d; n++)
474 boxv[d][n] += dt*t1[d][n];
476 /* We do NOT update the box vectors themselves here, since
477 * we need them for shifting later. It is instead done last
478 * in the update() routine.
481 /* Calculate the change relative to diagonal elements-
482 since it's perfectly ok for the off-diagonal ones to
483 be zero it doesn't make sense to check the change relative
487 change = fabs(dt*boxv[d][n]/box[d][d]);
489 if (change > maxchange)
496 if (maxchange > 0.01 && fplog)
500 "\nStep %s Warning: Pressure scaling more than 1%%. "
501 "This may mean your system\n is not yet equilibrated. "
502 "Use of Parrinello-Rahman pressure coupling during\n"
503 "equilibration can lead to simulation instability, "
504 "and is discouraged.\n",
505 gmx_step_str(step, buf));
509 preserve_box_shape(ir, box_rel, boxv);
511 mtmul(boxv, box, t1); /* t1=boxv * b' */
512 mmul(invbox, t1, t2);
513 mtmul(t2, invbox, M);
515 /* Determine the scaling matrix mu for the coordinates */
516 for (d = 0; d < DIM; d++)
518 for (n = 0; n <= d; n++)
520 t1[d][n] = box[d][n] + dt*boxv[d][n];
523 preserve_box_shape(ir, box_rel, t1);
524 /* t1 is the box at t+dt, determine mu as the relative change */
525 mmul_ur0(invbox, t1, mu);
528 void berendsen_pcoupl(FILE *fplog, gmx_int64_t step,
529 t_inputrec *ir, real dt, tensor pres, matrix box,
533 real scalar_pressure, xy_pressure, p_corr_z;
534 char *ptr, buf[STRLEN];
537 * Calculate the scaling matrix mu
541 for (d = 0; d < DIM; d++)
543 scalar_pressure += pres[d][d]/DIM;
546 xy_pressure += pres[d][d]/(DIM-1);
549 /* Pressure is now in bar, everywhere. */
550 #define factor(d, m) (ir->compress[d][m]*dt/ir->tau_p)
552 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
553 * necessary for triclinic scaling
559 for (d = 0; d < DIM; d++)
561 mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
564 case epctSEMIISOTROPIC:
565 for (d = 0; d < ZZ; d++)
567 mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
570 1.0 - factor(ZZ, ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
572 case epctANISOTROPIC:
573 for (d = 0; d < DIM; d++)
575 for (n = 0; n < DIM; n++)
577 mu[d][n] = (d == n ? 1.0 : 0.0)
578 -factor(d, n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
582 case epctSURFACETENSION:
583 /* ir->ref_p[0/1] is the reference surface-tension times *
584 * the number of surfaces */
585 if (ir->compress[ZZ][ZZ])
587 p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
591 /* when the compressibity is zero, set the pressure correction *
592 * in the z-direction to zero to get the correct surface tension */
595 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
596 for (d = 0; d < DIM-1; d++)
598 mu[d][d] = 1.0 + factor(d, d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
599 - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
603 gmx_fatal(FARGS, "Berendsen pressure coupling type %s not supported yet\n",
604 EPCOUPLTYPETYPE(ir->epct));
607 /* To fullfill the orientation restrictions on triclinic boxes
608 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
609 * the other elements of mu to first order.
611 mu[YY][XX] += mu[XX][YY];
612 mu[ZZ][XX] += mu[XX][ZZ];
613 mu[ZZ][YY] += mu[YY][ZZ];
620 pr_rvecs(debug, 0, "PC: pres ", pres, 3);
621 pr_rvecs(debug, 0, "PC: mu ", mu, 3);
624 if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 ||
625 mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01 ||
626 mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
629 sprintf(buf, "\nStep %s Warning: pressure scaling more than 1%%, "
631 gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
634 fprintf(fplog, "%s", buf);
636 fprintf(stderr, "%s", buf);
640 void berendsen_pscale(t_inputrec *ir, matrix mu,
641 matrix box, matrix box_rel,
642 int start, int nr_atoms,
643 rvec x[], unsigned short cFREEZE[],
646 ivec *nFreeze = ir->opts.nFreeze;
649 /* Scale the positions */
650 for (n = start; n < start+nr_atoms; n++)
659 x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
663 x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
667 x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
670 /* compute final boxlengths */
671 for (d = 0; d < DIM; d++)
673 box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
674 box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
675 box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
678 preserve_box_shape(ir, box_rel, box);
680 /* (un)shifting should NOT be done after this,
681 * since the box vectors might have changed
683 inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
686 void berendsen_tcoupl(t_inputrec *ir, gmx_ekindata_t *ekind, real dt)
690 real T, reft = 0, lll;
694 for (i = 0; (i < opts->ngtc); i++)
698 T = ekind->tcstat[i].T;
702 T = ekind->tcstat[i].Th;
705 if ((opts->tau_t[i] > 0) && (T > 0.0))
707 reft = max(0.0, opts->ref_t[i]);
708 lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
709 ekind->tcstat[i].lambda = max(min(lll, 1.25), 0.8);
713 ekind->tcstat[i].lambda = 1.0;
718 fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n",
719 i, T, ekind->tcstat[i].lambda);
724 void andersen_tcoupl(t_inputrec *ir, gmx_int64_t step,
725 const t_commrec *cr, const t_mdatoms *md, t_state *state, real rate, const gmx_bool *randomize, const real *boltzfac)
727 const int *gatindex = (DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
731 /* randomize the velocities of the selected particles */
733 for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */
735 int ng = gatindex ? gatindex[i] : i;
740 gc = md->cTC[i]; /* assign the atom to a temperature group if there are more than one */
744 if (ir->etc == etcANDERSEN)
752 gmx_rng_cycle_2uniform(step*2, ng, ir->andersen_seed, RND_SEED_ANDERSEN, uniform);
753 bRandomize = (uniform[0] < rate);
760 scal = sqrt(boltzfac[gc]*md->invmass[i]);
761 gmx_rng_cycle_3gaussian_table(step*2+1, ng, ir->andersen_seed, RND_SEED_ANDERSEN, gauss);
762 for (d = 0; d < DIM; d++)
764 state->v[i][d] = scal*gauss[d];
772 void nosehoover_tcoupl(t_grpopts *opts, gmx_ekindata_t *ekind, real dt,
773 double xi[], double vxi[], t_extmass *MassQ)
778 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
780 for (i = 0; (i < opts->ngtc); i++)
782 reft = max(0.0, opts->ref_t[i]);
784 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
785 xi[i] += dt*(oldvxi + vxi[i])*0.5;
789 t_state *init_bufstate(const t_state *template_state)
792 int nc = template_state->nhchainlength;
794 snew(state->nosehoover_xi, nc*template_state->ngtc);
795 snew(state->nosehoover_vxi, nc*template_state->ngtc);
796 snew(state->therm_integral, template_state->ngtc);
797 snew(state->nhpres_xi, nc*template_state->nnhpres);
798 snew(state->nhpres_vxi, nc*template_state->nnhpres);
803 void destroy_bufstate(t_state *state)
807 sfree(state->nosehoover_xi);
808 sfree(state->nosehoover_vxi);
809 sfree(state->therm_integral);
810 sfree(state->nhpres_xi);
811 sfree(state->nhpres_vxi);
815 void trotter_update(t_inputrec *ir, gmx_int64_t step, gmx_ekindata_t *ekind,
816 gmx_enerdata_t *enerd, t_state *state,
817 tensor vir, t_mdatoms *md,
818 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
821 int n, i, j, d, ntgrp, ngtc, gc = 0;
822 t_grp_tcstat *tcstat;
824 gmx_int64_t step_eff;
825 real ecorr, pcorr, dvdlcorr;
826 real bmass, qmass, reft, kT, dt, nd;
827 tensor dumpres, dumvir;
828 double *scalefac, dtc;
830 rvec sumv = {0, 0, 0}, consk;
833 if (trotter_seqno <= ettTSEQ2)
835 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
836 is actually the last half step from the previous step. Thus the first half step
837 actually corresponds to the n-1 step*/
845 bCouple = (ir->nsttcouple == 1 ||
846 do_per_step(step_eff+ir->nsttcouple, ir->nsttcouple));
848 trotter_seq = trotter_seqlist[trotter_seqno];
850 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
854 dtc = ir->nsttcouple*ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
855 opts = &(ir->opts); /* just for ease of referencing */
858 snew(scalefac, opts->ngtc);
859 for (i = 0; i < ngtc; i++)
863 /* execute the series of trotter updates specified in the trotterpart array */
865 for (i = 0; i < NTROTTERPARTS; i++)
867 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
868 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
877 switch (trotter_seq[i])
881 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir,
882 enerd->term[F_PDISPCORR], MassQ);
886 NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi,
887 state->nhpres_vxi, NULL, &(state->veta), MassQ, FALSE);
891 NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi,
892 state->nosehoover_vxi, scalefac, NULL, MassQ, (ir->eI == eiVV));
893 /* need to rescale the kinetic energies and velocities here. Could
894 scale the velocities later, but we need them scaled in order to
895 produce the correct outputs, so we'll scale them here. */
897 for (i = 0; i < ngtc; i++)
899 tcstat = &ekind->tcstat[i];
900 tcstat->vscale_nhc = scalefac[i];
901 tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
902 tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
904 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
905 /* but do we actually need the total? */
907 /* modify the velocities as well */
908 for (n = 0; n < md->homenr; n++)
910 if (md->cTC) /* does this conditional need to be here? is this always true?*/
914 for (d = 0; d < DIM; d++)
916 state->v[n][d] *= scalefac[gc];
921 for (d = 0; d < DIM; d++)
923 sumv[d] += (state->v[n][d])/md->invmass[n];
932 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
938 for (d = 0; d < DIM; d++)
940 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
942 fprintf(debug, "Conserved kappa: %15.8f %15.8f %15.8f\n", consk[0], consk[1], consk[2]);
950 extern void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit)
952 int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
953 t_grp_tcstat *tcstat;
955 real ecorr, pcorr, dvdlcorr;
956 real bmass, qmass, reft, kT, dt, ndj, nd;
957 tensor dumpres, dumvir;
959 opts = &(ir->opts); /* just for ease of referencing */
960 ngtc = ir->opts.ngtc;
961 nnhpres = state->nnhpres;
962 nh = state->nhchainlength;
968 snew(MassQ->Qinv, ngtc);
970 for (i = 0; (i < ngtc); i++)
972 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
974 MassQ->Qinv[i] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
978 MassQ->Qinv[i] = 0.0;
982 else if (EI_VV(ir->eI))
984 /* Set pressure variables */
988 if (state->vol0 == 0)
990 state->vol0 = det(state->box);
991 /* because we start by defining a fixed
992 compressibility, we need the volume at this
993 compressibility to solve the problem. */
997 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
998 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
999 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
1000 /* An alternate mass definition, from Tuckerman et al. */
1001 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1002 for (d = 0; d < DIM; d++)
1004 for (n = 0; n < DIM; n++)
1006 MassQ->Winvm[d][n] = PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
1007 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1008 before using MTTK for anisotropic states.*/
1011 /* Allocate space for thermostat variables */
1014 snew(MassQ->Qinv, ngtc*nh);
1017 /* now, set temperature variables */
1018 for (i = 0; i < ngtc; i++)
1020 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1022 reft = max(0.0, opts->ref_t[i]);
1025 for (j = 0; j < nh; j++)
1035 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
1041 for (j = 0; j < nh; j++)
1043 MassQ->Qinv[i*nh+j] = 0.0;
1050 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
1052 int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
1053 t_grp_tcstat *tcstat;
1055 real ecorr, pcorr, dvdlcorr;
1056 real bmass, qmass, reft, kT, dt, ndj, nd;
1057 tensor dumpres, dumvir;
1060 opts = &(ir->opts); /* just for ease of referencing */
1062 nnhpres = state->nnhpres;
1063 nh = state->nhchainlength;
1065 init_npt_masses(ir, state, MassQ, TRUE);
1067 /* first, initialize clear all the trotter calls */
1068 snew(trotter_seq, ettTSEQMAX);
1069 for (i = 0; i < ettTSEQMAX; i++)
1071 snew(trotter_seq[i], NTROTTERPARTS);
1072 for (j = 0; j < NTROTTERPARTS; j++)
1074 trotter_seq[i][j] = etrtNONE;
1076 trotter_seq[i][0] = etrtSKIPALL;
1081 /* no trotter calls, so we never use the values in the array.
1082 * We access them (so we need to define them, but ignore
1088 /* compute the kinetic energy by using the half step velocities or
1089 * the kinetic energies, depending on the order of the trotter calls */
1093 if (IR_NPT_TROTTER(ir))
1095 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1096 We start with the initial one. */
1097 /* first, a round that estimates veta. */
1098 trotter_seq[0][0] = etrtBAROV;
1100 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1102 /* The first half trotter update */
1103 trotter_seq[2][0] = etrtBAROV;
1104 trotter_seq[2][1] = etrtNHC;
1105 trotter_seq[2][2] = etrtBARONHC;
1107 /* The second half trotter update */
1108 trotter_seq[3][0] = etrtBARONHC;
1109 trotter_seq[3][1] = etrtNHC;
1110 trotter_seq[3][2] = etrtBAROV;
1112 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1115 else if (IR_NVT_TROTTER(ir))
1117 /* This is the easy version - there are only two calls, both the same.
1118 Otherwise, even easier -- no calls */
1119 trotter_seq[2][0] = etrtNHC;
1120 trotter_seq[3][0] = etrtNHC;
1122 else if (IR_NPH_TROTTER(ir))
1124 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1125 We start with the initial one. */
1126 /* first, a round that estimates veta. */
1127 trotter_seq[0][0] = etrtBAROV;
1129 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1131 /* The first half trotter update */
1132 trotter_seq[2][0] = etrtBAROV;
1133 trotter_seq[2][1] = etrtBARONHC;
1135 /* The second half trotter update */
1136 trotter_seq[3][0] = etrtBARONHC;
1137 trotter_seq[3][1] = etrtBAROV;
1139 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1142 else if (ir->eI == eiVVAK)
1144 if (IR_NPT_TROTTER(ir))
1146 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1147 We start with the initial one. */
1148 /* first, a round that estimates veta. */
1149 trotter_seq[0][0] = etrtBAROV;
1151 /* The first half trotter update, part 1 -- double update, because it commutes */
1152 trotter_seq[1][0] = etrtNHC;
1154 /* The first half trotter update, part 2 */
1155 trotter_seq[2][0] = etrtBAROV;
1156 trotter_seq[2][1] = etrtBARONHC;
1158 /* The second half trotter update, part 1 */
1159 trotter_seq[3][0] = etrtBARONHC;
1160 trotter_seq[3][1] = etrtBAROV;
1162 /* The second half trotter update */
1163 trotter_seq[4][0] = etrtNHC;
1165 else if (IR_NVT_TROTTER(ir))
1167 /* This is the easy version - there is only one call, both the same.
1168 Otherwise, even easier -- no calls */
1169 trotter_seq[1][0] = etrtNHC;
1170 trotter_seq[4][0] = etrtNHC;
1172 else if (IR_NPH_TROTTER(ir))
1174 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1175 We start with the initial one. */
1176 /* first, a round that estimates veta. */
1177 trotter_seq[0][0] = etrtBAROV;
1179 /* The first half trotter update, part 1 -- leave zero */
1180 trotter_seq[1][0] = etrtNHC;
1182 /* The first half trotter update, part 2 */
1183 trotter_seq[2][0] = etrtBAROV;
1184 trotter_seq[2][1] = etrtBARONHC;
1186 /* The second half trotter update, part 1 */
1187 trotter_seq[3][0] = etrtBARONHC;
1188 trotter_seq[3][1] = etrtBAROV;
1190 /* The second half trotter update -- blank for now */
1198 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
1201 snew(MassQ->QPinv, nnhpres*opts->nhchainlength);
1203 /* barostat temperature */
1204 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1206 reft = max(0.0, opts->ref_t[0]);
1208 for (i = 0; i < nnhpres; i++)
1210 for (j = 0; j < nh; j++)
1220 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1226 for (i = 0; i < nnhpres; i++)
1228 for (j = 0; j < nh; j++)
1230 MassQ->QPinv[i*nh+j] = 0.0;
1237 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1239 int i, j, nd, ndj, bmass, qmass, ngtcall;
1240 real ener_npt, reft, eta, kT, tau;
1243 real vol, dbaro, W, Q;
1244 int nh = state->nhchainlength;
1248 /* now we compute the contribution of the pressure to the conserved quantity*/
1250 if (ir->epc == epcMTTK)
1252 /* find the volume, and the kinetic energy of the volume */
1258 /* contribution from the pressure momenenta */
1259 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1261 /* contribution from the PV term */
1262 vol = det(state->box);
1263 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1266 case epctANISOTROPIC:
1270 case epctSURFACETENSION:
1273 case epctSEMIISOTROPIC:
1281 if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1283 /* add the energy from the barostat thermostat chain */
1284 for (i = 0; i < state->nnhpres; i++)
1287 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1288 ivxi = &state->nhpres_vxi[i*nh];
1289 ixi = &state->nhpres_xi[i*nh];
1290 iQinv = &(MassQ->QPinv[i*nh]);
1291 reft = max(ir->opts.ref_t[0], 0); /* using 'System' temperature */
1294 for (j = 0; j < nh; j++)
1298 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1299 /* contribution from the thermal variable of the NH chain */
1300 ener_npt += ixi[j]*kT;
1304 fprintf(debug, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i, j, ivxi[j], ixi[j]);
1312 for (i = 0; i < ir->opts.ngtc; i++)
1314 ixi = &state->nosehoover_xi[i*nh];
1315 ivxi = &state->nosehoover_vxi[i*nh];
1316 iQinv = &(MassQ->Qinv[i*nh]);
1318 nd = ir->opts.nrdf[i];
1319 reft = max(ir->opts.ref_t[i], 0);
1324 if (IR_NVT_TROTTER(ir))
1326 /* contribution from the thermal momenta of the NH chain */
1327 for (j = 0; j < nh; j++)
1331 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1332 /* contribution from the thermal variable of the NH chain */
1341 ener_npt += ndj*ixi[j]*kT;
1345 else /* Other non Trotter temperature NH control -- no chains yet. */
1347 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1348 ener_npt += nd*ixi[0]*kT;
1356 static real vrescale_gamdev(real ia,
1357 gmx_int64_t step, gmx_int64_t *count,
1358 gmx_int64_t seed1, gmx_int64_t seed2)
1359 /* Gamma distribution, adapted from numerical recipes */
1361 real am, e, s, v1, v2, x, y;
1372 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1374 v2 = 2.0*rnd[1] - 1.0;
1376 while (v1*v1 + v2*v2 > 1.0 ||
1377 v1*v1*GMX_REAL_MAX < 3.0*ia);
1378 /* The last check above ensures that both x (3.0 > 2.0 in s)
1379 * and the pre-factor for e do not go out of range.
1383 s = sqrt(2.0*am + 1.0);
1388 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1390 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1397 static real gaussian_count(gmx_int64_t step, gmx_int64_t *count,
1398 gmx_int64_t seed1, gmx_int64_t seed2)
1400 double rnd[2], x, y, r;
1404 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1405 x = 2.0*rnd[0] - 1.0;
1406 y = 2.0*rnd[1] - 1.0;
1409 while (r > 1.0 || r == 0.0);
1411 r = sqrt(-2.0*log(r)/r);
1416 static real vrescale_sumnoises(real nn,
1417 gmx_int64_t step, gmx_int64_t *count,
1418 gmx_int64_t seed1, gmx_int64_t seed2)
1421 * Returns the sum of nn independent gaussian noises squared
1422 * (i.e. equivalent to summing the square of the return values
1423 * of nn calls to gmx_rng_gaussian_real).
1425 const real ndeg_tol = 0.0001;
1428 if (nn < 2 + ndeg_tol)
1433 nn_int = (int)(nn + 0.5);
1435 if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1437 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);
1441 for (i = 0; i < nn_int; i++)
1443 gauss = gaussian_count(step, count, seed1, seed2);
1450 /* Use a gamma distribution for any real nn > 2 */
1451 r = 2.0*vrescale_gamdev(0.5*nn, step, count, seed1, seed2);
1457 static real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut,
1458 gmx_int64_t step, gmx_int64_t seed)
1461 * Generates a new value for the kinetic energy,
1462 * according to Bussi et al JCP (2007), Eq. (A7)
1463 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1464 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1465 * ndeg: number of degrees of freedom of the atoms to be thermalized
1466 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1468 /* rnd_count tracks the step-local state for the cycle random
1471 gmx_int64_t rnd_count = 0;
1472 real factor, rr, ekin_new;
1476 factor = exp(-1.0/taut);
1483 rr = gaussian_count(step, &rnd_count, seed, RND_SEED_VRESCALE);
1487 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, step, &rnd_count, seed, RND_SEED_VRESCALE) + rr*rr)/ndeg - kk) +
1488 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1493 void vrescale_tcoupl(t_inputrec *ir, gmx_int64_t step,
1494 gmx_ekindata_t *ekind, real dt,
1495 double therm_integral[])
1499 real Ek, Ek_ref1, Ek_ref, Ek_new;
1503 for (i = 0; (i < opts->ngtc); i++)
1507 Ek = trace(ekind->tcstat[i].ekinf);
1511 Ek = trace(ekind->tcstat[i].ekinh);
1514 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1516 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1517 Ek_ref = Ek_ref1*opts->nrdf[i];
1519 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i],
1523 /* Analytically Ek_new>=0, but we check for rounding errors */
1526 ekind->tcstat[i].lambda = 0.0;
1530 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1533 therm_integral[i] -= Ek_new - Ek;
1537 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1538 i, Ek_ref, Ek, Ek_new, ekind->tcstat[i].lambda);
1543 ekind->tcstat[i].lambda = 1.0;
1548 real vrescale_energy(t_grpopts *opts, double therm_integral[])
1554 for (i = 0; i < opts->ngtc; i++)
1556 ener += therm_integral[i];
1562 void rescale_velocities(gmx_ekindata_t *ekind, t_mdatoms *mdatoms,
1563 int start, int end, rvec v[])
1566 t_grp_tcstat *tcstat;
1567 unsigned short *cACC, *cTC;
1572 tcstat = ekind->tcstat;
1577 gstat = ekind->grpstat;
1578 cACC = mdatoms->cACC;
1582 for (n = start; n < end; n++)
1592 /* Only scale the velocity component relative to the COM velocity */
1593 rvec_sub(v[n], gstat[ga].u, vrel);
1594 lg = tcstat[gt].lambda;
1595 for (d = 0; d < DIM; d++)
1597 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1604 for (n = start; n < end; n++)
1610 lg = tcstat[gt].lambda;
1611 for (d = 0; d < DIM; d++)
1620 /* set target temperatures if we are annealing */
1621 void update_annealing_target_temp(t_grpopts *opts, real t)
1623 int i, j, n, npoints;
1624 real pert, thist = 0, x;
1626 for (i = 0; i < opts->ngtc; i++)
1628 npoints = opts->anneal_npoints[i];
1629 switch (opts->annealing[i])
1634 /* calculate time modulo the period */
1635 pert = opts->anneal_time[i][npoints-1];
1637 thist = t - n*pert; /* modulo time */
1638 /* Make sure rounding didn't get us outside the interval */
1639 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1648 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)", i, opts->ngtc, npoints);
1650 /* We are doing annealing for this group if we got here,
1651 * and we have the (relative) time as thist.
1652 * calculate target temp */
1654 while ((j < npoints-1) && (thist > (opts->anneal_time[i][j+1])))
1660 /* Found our position between points j and j+1.
1661 * Interpolate: x is the amount from j+1, (1-x) from point j
1662 * First treat possible jumps in temperature as a special case.
1664 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1666 opts->ref_t[i] = opts->anneal_temp[i][j+1];
1670 x = ((thist-opts->anneal_time[i][j])/
1671 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1672 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1677 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];