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.
42 #include "gromacs/legacyheaders/typedefs.h"
43 #include "gromacs/legacyheaders/types/commrec.h"
44 #include "gromacs/utility/smalloc.h"
45 #include "gromacs/legacyheaders/update.h"
46 #include "gromacs/math/vec.h"
47 #include "gromacs/legacyheaders/macros.h"
48 #include "gromacs/math/units.h"
49 #include "gromacs/legacyheaders/names.h"
50 #include "gromacs/utility/fatalerror.h"
51 #include "gromacs/legacyheaders/txtdump.h"
52 #include "gromacs/legacyheaders/nrnb.h"
53 #include "gromacs/random/random.h"
54 #include "gromacs/legacyheaders/update.h"
55 #include "gromacs/legacyheaders/mdrun.h"
57 #define NTROTTERPARTS 3
59 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
61 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
62 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
64 #define MAX_SUZUKI_YOSHIDA_NUM 5
65 #define SUZUKI_YOSHIDA_NUM 5
67 static const double sy_const_1[] = { 1. };
68 static const double sy_const_3[] = { 0.828981543588751, -0.657963087177502, 0.828981543588751 };
69 static const double sy_const_5[] = { 0.2967324292201065, 0.2967324292201065, -0.186929716880426, 0.2967324292201065, 0.2967324292201065 };
71 static const double* sy_const[] = {
81 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
85 {0.828981543588751,-0.657963087177502,0.828981543588751},
87 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
90 /* these integration routines are only referenced inside this file */
91 static void NHC_trotter(t_grpopts *opts, int nvar, gmx_ekindata_t *ekind, real dtfull,
92 double xi[], double vxi[], double scalefac[], real *veta, t_extmass *MassQ, gmx_bool bEkinAveVel)
95 /* general routine for both barostat and thermostat nose hoover chains */
97 int i, j, mi, mj, jmax;
98 double Ekin, Efac, reft, kT, nd;
100 t_grp_tcstat *tcstat;
105 int mstepsi, mstepsj;
106 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
107 int nh = opts->nhchainlength;
110 mstepsi = mstepsj = ns;
112 /* if scalefac is NULL, we are doing the NHC of the barostat */
115 if (scalefac == NULL)
120 for (i = 0; i < nvar; i++)
123 /* make it easier to iterate by selecting
124 out the sub-array that corresponds to this T group */
130 iQinv = &(MassQ->QPinv[i*nh]);
131 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
132 reft = max(0.0, opts->ref_t[0]);
133 Ekin = sqr(*veta)/MassQ->Winv;
137 iQinv = &(MassQ->Qinv[i*nh]);
138 tcstat = &ekind->tcstat[i];
140 reft = max(0.0, opts->ref_t[i]);
143 Ekin = 2*trace(tcstat->ekinf)*tcstat->ekinscalef_nhc;
147 Ekin = 2*trace(tcstat->ekinh)*tcstat->ekinscaleh_nhc;
152 for (mi = 0; mi < mstepsi; mi++)
154 for (mj = 0; mj < mstepsj; mj++)
156 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
157 dt = sy_const[ns][mj] * dtfull / mstepsi;
159 /* compute the thermal forces */
160 GQ[0] = iQinv[0]*(Ekin - nd*kT);
162 for (j = 0; j < nh-1; j++)
166 /* we actually don't need to update here if we save the
167 state of the GQ, but it's easier to just recompute*/
168 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
176 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
177 for (j = nh-1; j > 0; j--)
179 Efac = exp(-0.125*dt*ivxi[j]);
180 ivxi[j-1] = Efac*(ivxi[j-1]*Efac + 0.25*dt*GQ[j-1]);
183 Efac = exp(-0.5*dt*ivxi[0]);
194 /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
195 able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
196 think about this a bit more . . . */
198 GQ[0] = iQinv[0]*(Ekin - nd*kT);
200 /* update thermostat positions */
201 for (j = 0; j < nh; j++)
203 ixi[j] += 0.5*dt*ivxi[j];
206 for (j = 0; j < nh-1; j++)
208 Efac = exp(-0.125*dt*ivxi[j+1]);
209 ivxi[j] = Efac*(ivxi[j]*Efac + 0.25*dt*GQ[j]);
212 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
219 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
226 static void boxv_trotter(t_inputrec *ir, real *veta, real dt, tensor box,
227 gmx_ekindata_t *ekind, tensor vir, real pcorr, t_extmass *MassQ)
232 int i, j, d, n, nwall;
234 tensor Winvm, ekinmod, localpres;
236 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
237 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
240 if (ir->epct == epctSEMIISOTROPIC)
249 /* eta is in pure units. veta is in units of ps^-1. GW is in
250 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
251 taken to use only RATIOS of eta in updating the volume. */
253 /* we take the partial pressure tensors, modify the
254 kinetic energy tensor, and recovert to pressure */
256 if (ir->opts.nrdf[0] == 0)
258 gmx_fatal(FARGS, "Barostat is coupled to a T-group with no degrees of freedom\n");
260 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
261 alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
262 alpha *= ekind->tcstat[0].ekinscalef_nhc;
263 msmul(ekind->ekin, alpha, ekinmod);
264 /* for now, we use Elr = 0, because if you want to get it right, you
265 really should be using PME. Maybe print a warning? */
267 pscal = calc_pres(ir->ePBC, nwall, box, ekinmod, vir, localpres)+pcorr;
270 GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
276 * This file implements temperature and pressure coupling algorithms:
277 * For now only the Weak coupling and the modified weak coupling.
279 * Furthermore computation of pressure and temperature is done here
283 real calc_pres(int ePBC, int nwall, matrix box, tensor ekin, tensor vir,
289 if (ePBC == epbcNONE || (ePBC == epbcXY && nwall != 2))
295 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
296 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
300 fac = PRESFAC*2.0/det(box);
301 for (n = 0; (n < DIM); n++)
303 for (m = 0; (m < DIM); m++)
305 pres[n][m] = (ekin[n][m] - vir[n][m])*fac;
311 pr_rvecs(debug, 0, "PC: pres", pres, DIM);
312 pr_rvecs(debug, 0, "PC: ekin", ekin, DIM);
313 pr_rvecs(debug, 0, "PC: vir ", vir, DIM);
314 pr_rvecs(debug, 0, "PC: box ", box, DIM);
317 return trace(pres)/DIM;
320 real calc_temp(real ekin, real nrdf)
324 return (2.0*ekin)/(nrdf*BOLTZ);
332 void parrinellorahman_pcoupl(FILE *fplog, gmx_int64_t step,
333 t_inputrec *ir, real dt, tensor pres,
334 tensor box, tensor box_rel, tensor boxv,
335 tensor M, matrix mu, gmx_bool bFirstStep)
337 /* This doesn't do any coordinate updating. It just
338 * integrates the box vector equations from the calculated
339 * acceleration due to pressure difference. We also compute
340 * the tensor M which is used in update to couple the particle
341 * coordinates to the box vectors.
343 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
346 * M_nk = (h') * (h' * h + h' h) * h
348 * with the dots denoting time derivatives and h is the transformation from
349 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
350 * This also goes for the pressure and M tensors - they are transposed relative
351 * to ours. Our equation thus becomes:
354 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
356 * where b is the gromacs box matrix.
357 * Our box accelerations are given by
359 * b = vol/W inv(box') * (P-ref_P) (=h')
364 real vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
365 real atot, arel, change, maxchange, xy_pressure;
366 tensor invbox, pdiff, t1, t2;
370 m_inv_ur0(box, invbox);
374 /* Note that PRESFAC does not occur here.
375 * The pressure and compressibility always occur as a product,
376 * therefore the pressure unit drops out.
378 maxl = max(box[XX][XX], box[YY][YY]);
379 maxl = max(maxl, box[ZZ][ZZ]);
380 for (d = 0; d < DIM; d++)
382 for (n = 0; n < DIM; n++)
385 (4*M_PI*M_PI*ir->compress[d][n])/(3*ir->tau_p*ir->tau_p*maxl);
389 m_sub(pres, ir->ref_p, pdiff);
391 if (ir->epct == epctSURFACETENSION)
393 /* Unlike Berendsen coupling it might not be trivial to include a z
394 * pressure correction here? On the other hand we don't scale the
395 * box momentarily, but change accelerations, so it might not be crucial.
397 xy_pressure = 0.5*(pres[XX][XX]+pres[YY][YY]);
398 for (d = 0; d < ZZ; d++)
400 pdiff[d][d] = (xy_pressure-(pres[ZZ][ZZ]-ir->ref_p[d][d]/box[d][d]));
404 tmmul(invbox, pdiff, t1);
405 /* Move the off-diagonal elements of the 'force' to one side to ensure
406 * that we obey the box constraints.
408 for (d = 0; d < DIM; d++)
410 for (n = 0; n < d; n++)
412 t1[d][n] += t1[n][d];
419 case epctANISOTROPIC:
420 for (d = 0; d < DIM; d++)
422 for (n = 0; n <= d; n++)
424 t1[d][n] *= winv[d][n]*vol;
429 /* calculate total volume acceleration */
430 atot = box[XX][XX]*box[YY][YY]*t1[ZZ][ZZ]+
431 box[XX][XX]*t1[YY][YY]*box[ZZ][ZZ]+
432 t1[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
434 /* set all RELATIVE box accelerations equal, and maintain total V
436 for (d = 0; d < DIM; d++)
438 for (n = 0; n <= d; n++)
440 t1[d][n] = winv[0][0]*vol*arel*box[d][n];
444 case epctSEMIISOTROPIC:
445 case epctSURFACETENSION:
446 /* Note the correction to pdiff above for surftens. coupling */
448 /* calculate total XY volume acceleration */
449 atot = box[XX][XX]*t1[YY][YY]+t1[XX][XX]*box[YY][YY];
450 arel = atot/(2*box[XX][XX]*box[YY][YY]);
451 /* set RELATIVE XY box accelerations equal, and maintain total V
452 * change speed. Dont change the third box vector accelerations */
453 for (d = 0; d < ZZ; d++)
455 for (n = 0; n <= d; n++)
457 t1[d][n] = winv[d][n]*vol*arel*box[d][n];
460 for (n = 0; n < DIM; n++)
462 t1[ZZ][n] *= winv[d][n]*vol;
466 gmx_fatal(FARGS, "Parrinello-Rahman pressure coupling type %s "
467 "not supported yet\n", EPCOUPLTYPETYPE(ir->epct));
472 for (d = 0; d < DIM; d++)
474 for (n = 0; n <= d; n++)
476 boxv[d][n] += dt*t1[d][n];
478 /* We do NOT update the box vectors themselves here, since
479 * we need them for shifting later. It is instead done last
480 * in the update() routine.
483 /* Calculate the change relative to diagonal elements-
484 since it's perfectly ok for the off-diagonal ones to
485 be zero it doesn't make sense to check the change relative
489 change = fabs(dt*boxv[d][n]/box[d][d]);
491 if (change > maxchange)
498 if (maxchange > 0.01 && fplog)
502 "\nStep %s Warning: Pressure scaling more than 1%%. "
503 "This may mean your system\n is not yet equilibrated. "
504 "Use of Parrinello-Rahman pressure coupling during\n"
505 "equilibration can lead to simulation instability, "
506 "and is discouraged.\n",
507 gmx_step_str(step, buf));
511 preserve_box_shape(ir, box_rel, boxv);
513 mtmul(boxv, box, t1); /* t1=boxv * b' */
514 mmul(invbox, t1, t2);
515 mtmul(t2, invbox, M);
517 /* Determine the scaling matrix mu for the coordinates */
518 for (d = 0; d < DIM; d++)
520 for (n = 0; n <= d; n++)
522 t1[d][n] = box[d][n] + dt*boxv[d][n];
525 preserve_box_shape(ir, box_rel, t1);
526 /* t1 is the box at t+dt, determine mu as the relative change */
527 mmul_ur0(invbox, t1, mu);
530 void berendsen_pcoupl(FILE *fplog, gmx_int64_t step,
531 t_inputrec *ir, real dt, tensor pres, matrix box,
535 real scalar_pressure, xy_pressure, p_corr_z;
536 char *ptr, buf[STRLEN];
539 * Calculate the scaling matrix mu
543 for (d = 0; d < DIM; d++)
545 scalar_pressure += pres[d][d]/DIM;
548 xy_pressure += pres[d][d]/(DIM-1);
551 /* Pressure is now in bar, everywhere. */
552 #define factor(d, m) (ir->compress[d][m]*dt/ir->tau_p)
554 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
555 * necessary for triclinic scaling
561 for (d = 0; d < DIM; d++)
563 mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
566 case epctSEMIISOTROPIC:
567 for (d = 0; d < ZZ; d++)
569 mu[d][d] = 1.0 - factor(d, d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
572 1.0 - factor(ZZ, ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
574 case epctANISOTROPIC:
575 for (d = 0; d < DIM; d++)
577 for (n = 0; n < DIM; n++)
579 mu[d][n] = (d == n ? 1.0 : 0.0)
580 -factor(d, n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
584 case epctSURFACETENSION:
585 /* ir->ref_p[0/1] is the reference surface-tension times *
586 * the number of surfaces */
587 if (ir->compress[ZZ][ZZ])
589 p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
593 /* when the compressibity is zero, set the pressure correction *
594 * in the z-direction to zero to get the correct surface tension */
597 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
598 for (d = 0; d < DIM-1; d++)
600 mu[d][d] = 1.0 + factor(d, d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
601 - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
605 gmx_fatal(FARGS, "Berendsen pressure coupling type %s not supported yet\n",
606 EPCOUPLTYPETYPE(ir->epct));
609 /* To fullfill the orientation restrictions on triclinic boxes
610 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
611 * the other elements of mu to first order.
613 mu[YY][XX] += mu[XX][YY];
614 mu[ZZ][XX] += mu[XX][ZZ];
615 mu[ZZ][YY] += mu[YY][ZZ];
622 pr_rvecs(debug, 0, "PC: pres ", pres, 3);
623 pr_rvecs(debug, 0, "PC: mu ", mu, 3);
626 if (mu[XX][XX] < 0.99 || mu[XX][XX] > 1.01 ||
627 mu[YY][YY] < 0.99 || mu[YY][YY] > 1.01 ||
628 mu[ZZ][ZZ] < 0.99 || mu[ZZ][ZZ] > 1.01)
631 sprintf(buf, "\nStep %s Warning: pressure scaling more than 1%%, "
633 gmx_step_str(step, buf2), mu[XX][XX], mu[YY][YY], mu[ZZ][ZZ]);
636 fprintf(fplog, "%s", buf);
638 fprintf(stderr, "%s", buf);
642 void berendsen_pscale(t_inputrec *ir, matrix mu,
643 matrix box, matrix box_rel,
644 int start, int nr_atoms,
645 rvec x[], unsigned short cFREEZE[],
648 ivec *nFreeze = ir->opts.nFreeze;
651 /* Scale the positions */
652 for (n = start; n < start+nr_atoms; n++)
661 x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
665 x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
669 x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
672 /* compute final boxlengths */
673 for (d = 0; d < DIM; d++)
675 box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
676 box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
677 box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
680 preserve_box_shape(ir, box_rel, box);
682 /* (un)shifting should NOT be done after this,
683 * since the box vectors might have changed
685 inc_nrnb(nrnb, eNR_PCOUPL, nr_atoms);
688 void berendsen_tcoupl(t_inputrec *ir, gmx_ekindata_t *ekind, real dt)
692 real T, reft = 0, lll;
696 for (i = 0; (i < opts->ngtc); i++)
700 T = ekind->tcstat[i].T;
704 T = ekind->tcstat[i].Th;
707 if ((opts->tau_t[i] > 0) && (T > 0.0))
709 reft = max(0.0, opts->ref_t[i]);
710 lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
711 ekind->tcstat[i].lambda = max(min(lll, 1.25), 0.8);
715 ekind->tcstat[i].lambda = 1.0;
720 fprintf(debug, "TC: group %d: T: %g, Lambda: %g\n",
721 i, T, ekind->tcstat[i].lambda);
726 void andersen_tcoupl(t_inputrec *ir, gmx_int64_t step,
727 const t_commrec *cr, const t_mdatoms *md, t_state *state, real rate, const gmx_bool *randomize, const real *boltzfac)
729 const int *gatindex = (DOMAINDECOMP(cr) ? cr->dd->gatindex : NULL);
733 /* randomize the velocities of the selected particles */
735 for (i = 0; i < md->homenr; i++) /* now loop over the list of atoms */
737 int ng = gatindex ? gatindex[i] : i;
742 gc = md->cTC[i]; /* assign the atom to a temperature group if there are more than one */
746 if (ir->etc == etcANDERSEN)
754 gmx_rng_cycle_2uniform(step*2, ng, ir->andersen_seed, RND_SEED_ANDERSEN, uniform);
755 bRandomize = (uniform[0] < rate);
762 scal = sqrt(boltzfac[gc]*md->invmass[i]);
763 gmx_rng_cycle_3gaussian_table(step*2+1, ng, ir->andersen_seed, RND_SEED_ANDERSEN, gauss);
764 for (d = 0; d < DIM; d++)
766 state->v[i][d] = scal*gauss[d];
774 void nosehoover_tcoupl(t_grpopts *opts, gmx_ekindata_t *ekind, real dt,
775 double xi[], double vxi[], t_extmass *MassQ)
780 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
782 for (i = 0; (i < opts->ngtc); i++)
784 reft = max(0.0, opts->ref_t[i]);
786 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
787 xi[i] += dt*(oldvxi + vxi[i])*0.5;
791 t_state *init_bufstate(const t_state *template_state)
794 int nc = template_state->nhchainlength;
796 snew(state->nosehoover_xi, nc*template_state->ngtc);
797 snew(state->nosehoover_vxi, nc*template_state->ngtc);
798 snew(state->therm_integral, template_state->ngtc);
799 snew(state->nhpres_xi, nc*template_state->nnhpres);
800 snew(state->nhpres_vxi, nc*template_state->nnhpres);
805 void destroy_bufstate(t_state *state)
809 sfree(state->nosehoover_xi);
810 sfree(state->nosehoover_vxi);
811 sfree(state->therm_integral);
812 sfree(state->nhpres_xi);
813 sfree(state->nhpres_vxi);
817 void trotter_update(t_inputrec *ir, gmx_int64_t step, gmx_ekindata_t *ekind,
818 gmx_enerdata_t *enerd, t_state *state,
819 tensor vir, t_mdatoms *md,
820 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
823 int n, i, j, d, ntgrp, ngtc, gc = 0;
824 t_grp_tcstat *tcstat;
826 gmx_int64_t step_eff;
827 real ecorr, pcorr, dvdlcorr;
828 real bmass, qmass, reft, kT, dt, nd;
829 tensor dumpres, dumvir;
830 double *scalefac, dtc;
832 rvec sumv = {0, 0, 0}, consk;
835 if (trotter_seqno <= ettTSEQ2)
837 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
838 is actually the last half step from the previous step. Thus the first half step
839 actually corresponds to the n-1 step*/
847 bCouple = (ir->nsttcouple == 1 ||
848 do_per_step(step_eff+ir->nsttcouple, ir->nsttcouple));
850 trotter_seq = trotter_seqlist[trotter_seqno];
852 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
856 dtc = ir->nsttcouple*ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
857 opts = &(ir->opts); /* just for ease of referencing */
860 snew(scalefac, opts->ngtc);
861 for (i = 0; i < ngtc; i++)
865 /* execute the series of trotter updates specified in the trotterpart array */
867 for (i = 0; i < NTROTTERPARTS; i++)
869 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
870 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
879 switch (trotter_seq[i])
883 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir,
884 enerd->term[F_PDISPCORR], MassQ);
888 NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi,
889 state->nhpres_vxi, NULL, &(state->veta), MassQ, FALSE);
893 NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi,
894 state->nosehoover_vxi, scalefac, NULL, MassQ, (ir->eI == eiVV));
895 /* need to rescale the kinetic energies and velocities here. Could
896 scale the velocities later, but we need them scaled in order to
897 produce the correct outputs, so we'll scale them here. */
899 for (i = 0; i < ngtc; i++)
901 tcstat = &ekind->tcstat[i];
902 tcstat->vscale_nhc = scalefac[i];
903 tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
904 tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
906 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
907 /* but do we actually need the total? */
909 /* modify the velocities as well */
910 for (n = 0; n < md->homenr; n++)
912 if (md->cTC) /* does this conditional need to be here? is this always true?*/
916 for (d = 0; d < DIM; d++)
918 state->v[n][d] *= scalefac[gc];
923 for (d = 0; d < DIM; d++)
925 sumv[d] += (state->v[n][d])/md->invmass[n];
934 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
940 for (d = 0; d < DIM; d++)
942 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
944 fprintf(debug, "Conserved kappa: %15.8f %15.8f %15.8f\n", consk[0], consk[1], consk[2]);
952 extern void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit)
954 int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
955 t_grp_tcstat *tcstat;
957 real ecorr, pcorr, dvdlcorr;
958 real bmass, qmass, reft, kT, dt, ndj, nd;
959 tensor dumpres, dumvir;
961 opts = &(ir->opts); /* just for ease of referencing */
962 ngtc = ir->opts.ngtc;
963 nnhpres = state->nnhpres;
964 nh = state->nhchainlength;
970 snew(MassQ->Qinv, ngtc);
972 for (i = 0; (i < ngtc); i++)
974 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
976 MassQ->Qinv[i] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
980 MassQ->Qinv[i] = 0.0;
984 else if (EI_VV(ir->eI))
986 /* Set pressure variables */
990 if (state->vol0 == 0)
992 state->vol0 = det(state->box);
993 /* because we start by defining a fixed
994 compressibility, we need the volume at this
995 compressibility to solve the problem. */
999 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1000 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1001 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
1002 /* An alternate mass definition, from Tuckerman et al. */
1003 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1004 for (d = 0; d < DIM; d++)
1006 for (n = 0; n < DIM; n++)
1008 MassQ->Winvm[d][n] = PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
1009 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1010 before using MTTK for anisotropic states.*/
1013 /* Allocate space for thermostat variables */
1016 snew(MassQ->Qinv, ngtc*nh);
1019 /* now, set temperature variables */
1020 for (i = 0; i < ngtc; i++)
1022 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1024 reft = max(0.0, opts->ref_t[i]);
1027 for (j = 0; j < nh; j++)
1037 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
1043 for (j = 0; j < nh; j++)
1045 MassQ->Qinv[i*nh+j] = 0.0;
1052 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
1054 int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
1055 t_grp_tcstat *tcstat;
1057 real ecorr, pcorr, dvdlcorr;
1058 real bmass, qmass, reft, kT, dt, ndj, nd;
1059 tensor dumpres, dumvir;
1062 opts = &(ir->opts); /* just for ease of referencing */
1064 nnhpres = state->nnhpres;
1065 nh = state->nhchainlength;
1067 init_npt_masses(ir, state, MassQ, TRUE);
1069 /* first, initialize clear all the trotter calls */
1070 snew(trotter_seq, ettTSEQMAX);
1071 for (i = 0; i < ettTSEQMAX; i++)
1073 snew(trotter_seq[i], NTROTTERPARTS);
1074 for (j = 0; j < NTROTTERPARTS; j++)
1076 trotter_seq[i][j] = etrtNONE;
1078 trotter_seq[i][0] = etrtSKIPALL;
1083 /* no trotter calls, so we never use the values in the array.
1084 * We access them (so we need to define them, but ignore
1090 /* compute the kinetic energy by using the half step velocities or
1091 * the kinetic energies, depending on the order of the trotter calls */
1095 if (IR_NPT_TROTTER(ir))
1097 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1098 We start with the initial one. */
1099 /* first, a round that estimates veta. */
1100 trotter_seq[0][0] = etrtBAROV;
1102 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1104 /* The first half trotter update */
1105 trotter_seq[2][0] = etrtBAROV;
1106 trotter_seq[2][1] = etrtNHC;
1107 trotter_seq[2][2] = etrtBARONHC;
1109 /* The second half trotter update */
1110 trotter_seq[3][0] = etrtBARONHC;
1111 trotter_seq[3][1] = etrtNHC;
1112 trotter_seq[3][2] = etrtBAROV;
1114 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1117 else if (IR_NVT_TROTTER(ir))
1119 /* This is the easy version - there are only two calls, both the same.
1120 Otherwise, even easier -- no calls */
1121 trotter_seq[2][0] = etrtNHC;
1122 trotter_seq[3][0] = etrtNHC;
1124 else if (IR_NPH_TROTTER(ir))
1126 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1127 We start with the initial one. */
1128 /* first, a round that estimates veta. */
1129 trotter_seq[0][0] = etrtBAROV;
1131 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1133 /* The first half trotter update */
1134 trotter_seq[2][0] = etrtBAROV;
1135 trotter_seq[2][1] = etrtBARONHC;
1137 /* The second half trotter update */
1138 trotter_seq[3][0] = etrtBARONHC;
1139 trotter_seq[3][1] = etrtBAROV;
1141 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1144 else if (ir->eI == eiVVAK)
1146 if (IR_NPT_TROTTER(ir))
1148 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1149 We start with the initial one. */
1150 /* first, a round that estimates veta. */
1151 trotter_seq[0][0] = etrtBAROV;
1153 /* The first half trotter update, part 1 -- double update, because it commutes */
1154 trotter_seq[1][0] = etrtNHC;
1156 /* The first half trotter update, part 2 */
1157 trotter_seq[2][0] = etrtBAROV;
1158 trotter_seq[2][1] = etrtBARONHC;
1160 /* The second half trotter update, part 1 */
1161 trotter_seq[3][0] = etrtBARONHC;
1162 trotter_seq[3][1] = etrtBAROV;
1164 /* The second half trotter update */
1165 trotter_seq[4][0] = etrtNHC;
1167 else if (IR_NVT_TROTTER(ir))
1169 /* This is the easy version - there is only one call, both the same.
1170 Otherwise, even easier -- no calls */
1171 trotter_seq[1][0] = etrtNHC;
1172 trotter_seq[4][0] = etrtNHC;
1174 else if (IR_NPH_TROTTER(ir))
1176 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1177 We start with the initial one. */
1178 /* first, a round that estimates veta. */
1179 trotter_seq[0][0] = etrtBAROV;
1181 /* The first half trotter update, part 1 -- leave zero */
1182 trotter_seq[1][0] = etrtNHC;
1184 /* The first half trotter update, part 2 */
1185 trotter_seq[2][0] = etrtBAROV;
1186 trotter_seq[2][1] = etrtBARONHC;
1188 /* The second half trotter update, part 1 */
1189 trotter_seq[3][0] = etrtBARONHC;
1190 trotter_seq[3][1] = etrtBAROV;
1192 /* The second half trotter update -- blank for now */
1200 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
1203 snew(MassQ->QPinv, nnhpres*opts->nhchainlength);
1205 /* barostat temperature */
1206 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1208 reft = max(0.0, opts->ref_t[0]);
1210 for (i = 0; i < nnhpres; i++)
1212 for (j = 0; j < nh; j++)
1222 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1228 for (i = 0; i < nnhpres; i++)
1230 for (j = 0; j < nh; j++)
1232 MassQ->QPinv[i*nh+j] = 0.0;
1239 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1241 int i, j, nd, ndj, bmass, qmass, ngtcall;
1242 real ener_npt, reft, eta, kT, tau;
1245 real vol, dbaro, W, Q;
1246 int nh = state->nhchainlength;
1250 /* now we compute the contribution of the pressure to the conserved quantity*/
1252 if (ir->epc == epcMTTK)
1254 /* find the volume, and the kinetic energy of the volume */
1260 /* contribution from the pressure momenenta */
1261 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1263 /* contribution from the PV term */
1264 vol = det(state->box);
1265 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1268 case epctANISOTROPIC:
1272 case epctSURFACETENSION:
1275 case epctSEMIISOTROPIC:
1283 if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1285 /* add the energy from the barostat thermostat chain */
1286 for (i = 0; i < state->nnhpres; i++)
1289 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1290 ivxi = &state->nhpres_vxi[i*nh];
1291 ixi = &state->nhpres_xi[i*nh];
1292 iQinv = &(MassQ->QPinv[i*nh]);
1293 reft = max(ir->opts.ref_t[0], 0); /* using 'System' temperature */
1296 for (j = 0; j < nh; j++)
1300 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1301 /* contribution from the thermal variable of the NH chain */
1302 ener_npt += ixi[j]*kT;
1306 fprintf(debug, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i, j, ivxi[j], ixi[j]);
1314 for (i = 0; i < ir->opts.ngtc; i++)
1316 ixi = &state->nosehoover_xi[i*nh];
1317 ivxi = &state->nosehoover_vxi[i*nh];
1318 iQinv = &(MassQ->Qinv[i*nh]);
1320 nd = ir->opts.nrdf[i];
1321 reft = max(ir->opts.ref_t[i], 0);
1326 if (IR_NVT_TROTTER(ir))
1328 /* contribution from the thermal momenta of the NH chain */
1329 for (j = 0; j < nh; j++)
1333 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1334 /* contribution from the thermal variable of the NH chain */
1343 ener_npt += ndj*ixi[j]*kT;
1347 else /* Other non Trotter temperature NH control -- no chains yet. */
1349 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1350 ener_npt += nd*ixi[0]*kT;
1358 static real vrescale_gamdev(real ia,
1359 gmx_int64_t step, gmx_int64_t *count,
1360 gmx_int64_t seed1, gmx_int64_t seed2)
1361 /* Gamma distribution, adapted from numerical recipes */
1363 real am, e, s, v1, v2, x, y;
1374 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1376 v2 = 2.0*rnd[1] - 1.0;
1378 while (v1*v1 + v2*v2 > 1.0 ||
1379 v1*v1*GMX_REAL_MAX < 3.0*ia);
1380 /* The last check above ensures that both x (3.0 > 2.0 in s)
1381 * and the pre-factor for e do not go out of range.
1385 s = sqrt(2.0*am + 1.0);
1390 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1392 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1399 static real gaussian_count(gmx_int64_t step, gmx_int64_t *count,
1400 gmx_int64_t seed1, gmx_int64_t seed2)
1402 double rnd[2], x, y, r;
1406 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1407 x = 2.0*rnd[0] - 1.0;
1408 y = 2.0*rnd[1] - 1.0;
1411 while (r > 1.0 || r == 0.0);
1413 r = sqrt(-2.0*log(r)/r);
1418 static real vrescale_sumnoises(real nn,
1419 gmx_int64_t step, gmx_int64_t *count,
1420 gmx_int64_t seed1, gmx_int64_t seed2)
1423 * Returns the sum of nn independent gaussian noises squared
1424 * (i.e. equivalent to summing the square of the return values
1425 * of nn calls to gmx_rng_gaussian_real).
1427 const real ndeg_tol = 0.0001;
1430 if (nn < 2 + ndeg_tol)
1435 nn_int = (int)(nn + 0.5);
1437 if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1439 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);
1443 for (i = 0; i < nn_int; i++)
1445 gauss = gaussian_count(step, count, seed1, seed2);
1452 /* Use a gamma distribution for any real nn > 2 */
1453 r = 2.0*vrescale_gamdev(0.5*nn, step, count, seed1, seed2);
1459 static real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut,
1460 gmx_int64_t step, gmx_int64_t seed)
1463 * Generates a new value for the kinetic energy,
1464 * according to Bussi et al JCP (2007), Eq. (A7)
1465 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1466 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1467 * ndeg: number of degrees of freedom of the atoms to be thermalized
1468 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1470 /* rnd_count tracks the step-local state for the cycle random
1473 gmx_int64_t rnd_count = 0;
1474 real factor, rr, ekin_new;
1478 factor = exp(-1.0/taut);
1485 rr = gaussian_count(step, &rnd_count, seed, RND_SEED_VRESCALE);
1489 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, step, &rnd_count, seed, RND_SEED_VRESCALE) + rr*rr)/ndeg - kk) +
1490 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1495 void vrescale_tcoupl(t_inputrec *ir, gmx_int64_t step,
1496 gmx_ekindata_t *ekind, real dt,
1497 double therm_integral[])
1501 real Ek, Ek_ref1, Ek_ref, Ek_new;
1505 for (i = 0; (i < opts->ngtc); i++)
1509 Ek = trace(ekind->tcstat[i].ekinf);
1513 Ek = trace(ekind->tcstat[i].ekinh);
1516 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1518 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1519 Ek_ref = Ek_ref1*opts->nrdf[i];
1521 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i],
1525 /* Analytically Ek_new>=0, but we check for rounding errors */
1528 ekind->tcstat[i].lambda = 0.0;
1532 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1535 therm_integral[i] -= Ek_new - Ek;
1539 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1540 i, Ek_ref, Ek, Ek_new, ekind->tcstat[i].lambda);
1545 ekind->tcstat[i].lambda = 1.0;
1550 real vrescale_energy(t_grpopts *opts, double therm_integral[])
1556 for (i = 0; i < opts->ngtc; i++)
1558 ener += therm_integral[i];
1564 void rescale_velocities(gmx_ekindata_t *ekind, t_mdatoms *mdatoms,
1565 int start, int end, rvec v[])
1568 t_grp_tcstat *tcstat;
1569 unsigned short *cACC, *cTC;
1574 tcstat = ekind->tcstat;
1579 gstat = ekind->grpstat;
1580 cACC = mdatoms->cACC;
1584 for (n = start; n < end; n++)
1594 /* Only scale the velocity component relative to the COM velocity */
1595 rvec_sub(v[n], gstat[ga].u, vrel);
1596 lg = tcstat[gt].lambda;
1597 for (d = 0; d < DIM; d++)
1599 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1606 for (n = start; n < end; n++)
1612 lg = tcstat[gt].lambda;
1613 for (d = 0; d < DIM; d++)
1622 /* set target temperatures if we are annealing */
1623 void update_annealing_target_temp(t_grpopts *opts, real t)
1625 int i, j, n, npoints;
1626 real pert, thist = 0, x;
1628 for (i = 0; i < opts->ngtc; i++)
1630 npoints = opts->anneal_npoints[i];
1631 switch (opts->annealing[i])
1636 /* calculate time modulo the period */
1637 pert = opts->anneal_time[i][npoints-1];
1639 thist = t - n*pert; /* modulo time */
1640 /* Make sure rounding didn't get us outside the interval */
1641 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1650 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)", i, opts->ngtc, npoints);
1652 /* We are doing annealing for this group if we got here,
1653 * and we have the (relative) time as thist.
1654 * calculate target temp */
1656 while ((j < npoints-1) && (thist > (opts->anneal_time[i][j+1])))
1662 /* Found our position between points j and j+1.
1663 * Interpolate: x is the amount from j+1, (1-x) from point j
1664 * First treat possible jumps in temperature as a special case.
1666 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1668 opts->ref_t[i] = opts->anneal_temp[i][j+1];
1672 x = ((thist-opts->anneal_time[i][j])/
1673 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1674 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1679 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];