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.
43 #include "gromacs/legacyheaders/macros.h"
44 #include "gromacs/legacyheaders/mdrun.h"
45 #include "gromacs/legacyheaders/names.h"
46 #include "gromacs/legacyheaders/nrnb.h"
47 #include "gromacs/legacyheaders/txtdump.h"
48 #include "gromacs/legacyheaders/typedefs.h"
49 #include "gromacs/legacyheaders/update.h"
50 #include "gromacs/legacyheaders/types/commrec.h"
51 #include "gromacs/math/units.h"
52 #include "gromacs/math/vec.h"
53 #include "gromacs/random/random.h"
54 #include "gromacs/utility/fatalerror.h"
55 #include "gromacs/utility/smalloc.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 */
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 = std::max<real>(0, opts->ref_t[0]);
133 Ekin = sqr(*veta)/MassQ->Winv;
137 iQinv = &(MassQ->Qinv[i*nh]);
138 tcstat = &ekind->tcstat[i];
140 reft = std::max<real>(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)
234 tensor 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 = std::max(box[XX][XX], box[YY][YY]);
379 maxl = std::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;
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 = std::max<real>(0, opts->ref_t[i]);
710 lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
711 ekind->tcstat[i].lambda = std::max<real>(std::min<real>(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 = std::max<real>(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, d, ngtc, gc = 0;
824 t_grp_tcstat *tcstat;
826 gmx_int64_t step_eff;
828 double *scalefac, dtc;
830 rvec sumv = {0, 0, 0};
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, ngtc, nh;
954 real reft, kT, ndj, nd;
956 opts = &(ir->opts); /* just for ease of referencing */
957 ngtc = ir->opts.ngtc;
958 nh = state->nhchainlength;
964 snew(MassQ->Qinv, ngtc);
966 for (i = 0; (i < ngtc); i++)
968 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
970 MassQ->Qinv[i] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
974 MassQ->Qinv[i] = 0.0;
978 else if (EI_VV(ir->eI))
980 /* Set pressure variables */
984 if (state->vol0 == 0)
986 state->vol0 = det(state->box);
987 /* because we start by defining a fixed
988 compressibility, we need the volume at this
989 compressibility to solve the problem. */
993 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
994 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
995 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
996 /* An alternate mass definition, from Tuckerman et al. */
997 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
998 for (d = 0; d < DIM; d++)
1000 for (n = 0; n < DIM; n++)
1002 MassQ->Winvm[d][n] = PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
1003 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1004 before using MTTK for anisotropic states.*/
1007 /* Allocate space for thermostat variables */
1010 snew(MassQ->Qinv, ngtc*nh);
1013 /* now, set temperature variables */
1014 for (i = 0; i < ngtc; i++)
1016 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1018 reft = std::max<real>(0, opts->ref_t[i]);
1021 for (j = 0; j < nh; j++)
1031 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
1036 for (j = 0; j < nh; j++)
1038 MassQ->Qinv[i*nh+j] = 0.0;
1045 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
1047 int i, j, nnhpres, nh;
1049 real bmass, qmass, reft, kT;
1052 opts = &(ir->opts); /* just for ease of referencing */
1053 nnhpres = state->nnhpres;
1054 nh = state->nhchainlength;
1056 init_npt_masses(ir, state, MassQ, TRUE);
1058 /* first, initialize clear all the trotter calls */
1059 snew(trotter_seq, ettTSEQMAX);
1060 for (i = 0; i < ettTSEQMAX; i++)
1062 snew(trotter_seq[i], NTROTTERPARTS);
1063 for (j = 0; j < NTROTTERPARTS; j++)
1065 trotter_seq[i][j] = etrtNONE;
1067 trotter_seq[i][0] = etrtSKIPALL;
1072 /* no trotter calls, so we never use the values in the array.
1073 * We access them (so we need to define them, but ignore
1079 /* compute the kinetic energy by using the half step velocities or
1080 * the kinetic energies, depending on the order of the trotter calls */
1084 if (IR_NPT_TROTTER(ir))
1086 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1087 We start with the initial one. */
1088 /* first, a round that estimates veta. */
1089 trotter_seq[0][0] = etrtBAROV;
1091 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1093 /* The first half trotter update */
1094 trotter_seq[2][0] = etrtBAROV;
1095 trotter_seq[2][1] = etrtNHC;
1096 trotter_seq[2][2] = etrtBARONHC;
1098 /* The second half trotter update */
1099 trotter_seq[3][0] = etrtBARONHC;
1100 trotter_seq[3][1] = etrtNHC;
1101 trotter_seq[3][2] = etrtBAROV;
1103 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1106 else if (IR_NVT_TROTTER(ir))
1108 /* This is the easy version - there are only two calls, both the same.
1109 Otherwise, even easier -- no calls */
1110 trotter_seq[2][0] = etrtNHC;
1111 trotter_seq[3][0] = etrtNHC;
1113 else if (IR_NPH_TROTTER(ir))
1115 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1116 We start with the initial one. */
1117 /* first, a round that estimates veta. */
1118 trotter_seq[0][0] = etrtBAROV;
1120 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1122 /* The first half trotter update */
1123 trotter_seq[2][0] = etrtBAROV;
1124 trotter_seq[2][1] = etrtBARONHC;
1126 /* The second half trotter update */
1127 trotter_seq[3][0] = etrtBARONHC;
1128 trotter_seq[3][1] = etrtBAROV;
1130 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1133 else if (ir->eI == eiVVAK)
1135 if (IR_NPT_TROTTER(ir))
1137 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1138 We start with the initial one. */
1139 /* first, a round that estimates veta. */
1140 trotter_seq[0][0] = etrtBAROV;
1142 /* The first half trotter update, part 1 -- double update, because it commutes */
1143 trotter_seq[1][0] = etrtNHC;
1145 /* The first half trotter update, part 2 */
1146 trotter_seq[2][0] = etrtBAROV;
1147 trotter_seq[2][1] = etrtBARONHC;
1149 /* The second half trotter update, part 1 */
1150 trotter_seq[3][0] = etrtBARONHC;
1151 trotter_seq[3][1] = etrtBAROV;
1153 /* The second half trotter update */
1154 trotter_seq[4][0] = etrtNHC;
1156 else if (IR_NVT_TROTTER(ir))
1158 /* This is the easy version - there is only one call, both the same.
1159 Otherwise, even easier -- no calls */
1160 trotter_seq[1][0] = etrtNHC;
1161 trotter_seq[4][0] = etrtNHC;
1163 else if (IR_NPH_TROTTER(ir))
1165 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1166 We start with the initial one. */
1167 /* first, a round that estimates veta. */
1168 trotter_seq[0][0] = etrtBAROV;
1170 /* The first half trotter update, part 1 -- leave zero */
1171 trotter_seq[1][0] = etrtNHC;
1173 /* The first half trotter update, part 2 */
1174 trotter_seq[2][0] = etrtBAROV;
1175 trotter_seq[2][1] = etrtBARONHC;
1177 /* The second half trotter update, part 1 */
1178 trotter_seq[3][0] = etrtBARONHC;
1179 trotter_seq[3][1] = etrtBAROV;
1181 /* The second half trotter update -- blank for now */
1189 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
1192 snew(MassQ->QPinv, nnhpres*opts->nhchainlength);
1194 /* barostat temperature */
1195 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1197 reft = std::max<real>(0, opts->ref_t[0]);
1199 for (i = 0; i < nnhpres; i++)
1201 for (j = 0; j < nh; j++)
1211 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1217 for (i = 0; i < nnhpres; i++)
1219 for (j = 0; j < nh; j++)
1221 MassQ->QPinv[i*nh+j] = 0.0;
1228 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1232 real ener_npt, reft, kT;
1236 int nh = state->nhchainlength;
1240 /* now we compute the contribution of the pressure to the conserved quantity*/
1242 if (ir->epc == epcMTTK)
1244 /* find the volume, and the kinetic energy of the volume */
1250 /* contribution from the pressure momenenta */
1251 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1253 /* contribution from the PV term */
1254 vol = det(state->box);
1255 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1258 case epctANISOTROPIC:
1262 case epctSURFACETENSION:
1265 case epctSEMIISOTROPIC:
1273 if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1275 /* add the energy from the barostat thermostat chain */
1276 for (i = 0; i < state->nnhpres; i++)
1279 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1280 ivxi = &state->nhpres_vxi[i*nh];
1281 ixi = &state->nhpres_xi[i*nh];
1282 iQinv = &(MassQ->QPinv[i*nh]);
1283 reft = std::max<real>(ir->opts.ref_t[0], 0.0); /* using 'System' temperature */
1286 for (j = 0; j < nh; j++)
1290 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1291 /* contribution from the thermal variable of the NH chain */
1292 ener_npt += ixi[j]*kT;
1296 fprintf(debug, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i, j, ivxi[j], ixi[j]);
1304 for (i = 0; i < ir->opts.ngtc; i++)
1306 ixi = &state->nosehoover_xi[i*nh];
1307 ivxi = &state->nosehoover_vxi[i*nh];
1308 iQinv = &(MassQ->Qinv[i*nh]);
1310 nd = ir->opts.nrdf[i];
1311 reft = std::max<real>(ir->opts.ref_t[i], 0);
1316 if (IR_NVT_TROTTER(ir))
1318 /* contribution from the thermal momenta of the NH chain */
1319 for (j = 0; j < nh; j++)
1323 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1324 /* contribution from the thermal variable of the NH chain */
1333 ener_npt += ndj*ixi[j]*kT;
1337 else /* Other non Trotter temperature NH control -- no chains yet. */
1339 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1340 ener_npt += nd*ixi[0]*kT;
1348 static real vrescale_gamdev(real ia,
1349 gmx_int64_t step, gmx_int64_t *count,
1350 gmx_int64_t seed1, gmx_int64_t seed2)
1351 /* Gamma distribution, adapted from numerical recipes */
1353 real am, e, s, v1, v2, x, y;
1364 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1366 v2 = 2.0*rnd[1] - 1.0;
1368 while (v1*v1 + v2*v2 > 1.0 ||
1369 v1*v1*GMX_REAL_MAX < 3.0*ia);
1370 /* The last check above ensures that both x (3.0 > 2.0 in s)
1371 * and the pre-factor for e do not go out of range.
1375 s = sqrt(2.0*am + 1.0);
1380 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1382 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1389 static real gaussian_count(gmx_int64_t step, gmx_int64_t *count,
1390 gmx_int64_t seed1, gmx_int64_t seed2)
1392 double rnd[2], x, y, r;
1396 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1397 x = 2.0*rnd[0] - 1.0;
1398 y = 2.0*rnd[1] - 1.0;
1401 while (r > 1.0 || r == 0.0);
1403 r = sqrt(-2.0*log(r)/r);
1408 static real vrescale_sumnoises(real nn,
1409 gmx_int64_t step, gmx_int64_t *count,
1410 gmx_int64_t seed1, gmx_int64_t seed2)
1413 * Returns the sum of nn independent gaussian noises squared
1414 * (i.e. equivalent to summing the square of the return values
1415 * of nn calls to gmx_rng_gaussian_real).
1417 const real ndeg_tol = 0.0001;
1420 if (nn < 2 + ndeg_tol)
1425 nn_int = (int)(nn + 0.5);
1427 if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1429 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);
1433 for (i = 0; i < nn_int; i++)
1435 gauss = gaussian_count(step, count, seed1, seed2);
1442 /* Use a gamma distribution for any real nn > 2 */
1443 r = 2.0*vrescale_gamdev(0.5*nn, step, count, seed1, seed2);
1449 static real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut,
1450 gmx_int64_t step, gmx_int64_t seed)
1453 * Generates a new value for the kinetic energy,
1454 * according to Bussi et al JCP (2007), Eq. (A7)
1455 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1456 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1457 * ndeg: number of degrees of freedom of the atoms to be thermalized
1458 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1460 /* rnd_count tracks the step-local state for the cycle random
1463 gmx_int64_t rnd_count = 0;
1464 real factor, rr, ekin_new;
1468 factor = exp(-1.0/taut);
1475 rr = gaussian_count(step, &rnd_count, seed, RND_SEED_VRESCALE);
1479 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, step, &rnd_count, seed, RND_SEED_VRESCALE) + rr*rr)/ndeg - kk) +
1480 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1485 void vrescale_tcoupl(t_inputrec *ir, gmx_int64_t step,
1486 gmx_ekindata_t *ekind, real dt,
1487 double therm_integral[])
1491 real Ek, Ek_ref1, Ek_ref, Ek_new;
1495 for (i = 0; (i < opts->ngtc); i++)
1499 Ek = trace(ekind->tcstat[i].ekinf);
1503 Ek = trace(ekind->tcstat[i].ekinh);
1506 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1508 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1509 Ek_ref = Ek_ref1*opts->nrdf[i];
1511 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i],
1515 /* Analytically Ek_new>=0, but we check for rounding errors */
1518 ekind->tcstat[i].lambda = 0.0;
1522 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1525 therm_integral[i] -= Ek_new - Ek;
1529 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1530 i, Ek_ref, Ek, Ek_new, ekind->tcstat[i].lambda);
1535 ekind->tcstat[i].lambda = 1.0;
1540 real vrescale_energy(t_grpopts *opts, double therm_integral[])
1546 for (i = 0; i < opts->ngtc; i++)
1548 ener += therm_integral[i];
1554 void rescale_velocities(gmx_ekindata_t *ekind, t_mdatoms *mdatoms,
1555 int start, int end, rvec v[])
1558 t_grp_tcstat *tcstat;
1559 unsigned short *cACC, *cTC;
1564 tcstat = ekind->tcstat;
1569 gstat = ekind->grpstat;
1570 cACC = mdatoms->cACC;
1574 for (n = start; n < end; n++)
1584 /* Only scale the velocity component relative to the COM velocity */
1585 rvec_sub(v[n], gstat[ga].u, vrel);
1586 lg = tcstat[gt].lambda;
1587 for (d = 0; d < DIM; d++)
1589 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1596 for (n = start; n < end; n++)
1602 lg = tcstat[gt].lambda;
1603 for (d = 0; d < DIM; d++)
1612 /* set target temperatures if we are annealing */
1613 void update_annealing_target_temp(t_grpopts *opts, real t)
1615 int i, j, n, npoints;
1616 real pert, thist = 0, x;
1618 for (i = 0; i < opts->ngtc; i++)
1620 npoints = opts->anneal_npoints[i];
1621 switch (opts->annealing[i])
1626 /* calculate time modulo the period */
1627 pert = opts->anneal_time[i][npoints-1];
1628 n = static_cast<int>(t / pert);
1629 thist = t - n*pert; /* modulo time */
1630 /* Make sure rounding didn't get us outside the interval */
1631 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1640 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)", i, opts->ngtc, npoints);
1642 /* We are doing annealing for this group if we got here,
1643 * and we have the (relative) time as thist.
1644 * calculate target temp */
1646 while ((j < npoints-1) && (thist > (opts->anneal_time[i][j+1])))
1652 /* Found our position between points j and j+1.
1653 * Interpolate: x is the amount from j+1, (1-x) from point j
1654 * First treat possible jumps in temperature as a special case.
1656 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1658 opts->ref_t[i] = opts->anneal_temp[i][j+1];
1662 x = ((thist-opts->anneal_time[i][j])/
1663 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1664 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1669 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];