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,2015, 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 "types/commrec.h"
44 #include "gromacs/utility/smalloc.h"
50 #include "gmx_fatal.h"
53 #include "gromacs/random/random.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 == etcANDERSENMASSIVE)
748 /* Randomize particle always */
753 /* Randomize particle probabilistically */
756 gmx_rng_cycle_2uniform(step*2, ng, ir->andersen_seed, RND_SEED_ANDERSEN, uniform);
757 bRandomize = (uniform[0] < rate);
764 scal = sqrt(boltzfac[gc]*md->invmass[i]);
765 gmx_rng_cycle_3gaussian_table(step*2+1, ng, ir->andersen_seed, RND_SEED_ANDERSEN, gauss);
766 for (d = 0; d < DIM; d++)
768 state->v[i][d] = scal*gauss[d];
776 void nosehoover_tcoupl(t_grpopts *opts, gmx_ekindata_t *ekind, real dt,
777 double xi[], double vxi[], t_extmass *MassQ)
782 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
784 for (i = 0; (i < opts->ngtc); i++)
786 reft = max(0.0, opts->ref_t[i]);
788 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
789 xi[i] += dt*(oldvxi + vxi[i])*0.5;
793 t_state *init_bufstate(const t_state *template_state)
796 int nc = template_state->nhchainlength;
798 snew(state->nosehoover_xi, nc*template_state->ngtc);
799 snew(state->nosehoover_vxi, nc*template_state->ngtc);
800 snew(state->therm_integral, template_state->ngtc);
801 snew(state->nhpres_xi, nc*template_state->nnhpres);
802 snew(state->nhpres_vxi, nc*template_state->nnhpres);
807 void destroy_bufstate(t_state *state)
811 sfree(state->nosehoover_xi);
812 sfree(state->nosehoover_vxi);
813 sfree(state->therm_integral);
814 sfree(state->nhpres_xi);
815 sfree(state->nhpres_vxi);
819 void trotter_update(t_inputrec *ir, gmx_int64_t step, gmx_ekindata_t *ekind,
820 gmx_enerdata_t *enerd, t_state *state,
821 tensor vir, t_mdatoms *md,
822 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
825 int n, i, j, d, ntgrp, ngtc, gc = 0, t;
826 t_grp_tcstat *tcstat;
828 gmx_int64_t step_eff;
829 real ecorr, pcorr, dvdlcorr;
830 real bmass, qmass, reft, kT, dt, nd;
831 tensor dumpres, dumvir;
832 double *scalefac, dtc;
834 rvec sumv = {0, 0, 0}, consk;
837 if (trotter_seqno <= ettTSEQ2)
839 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
840 is actually the last half step from the previous step. Thus the first half step
841 actually corresponds to the n-1 step*/
849 bCouple = (ir->nsttcouple == 1 ||
850 do_per_step(step_eff+ir->nsttcouple, ir->nsttcouple));
852 trotter_seq = trotter_seqlist[trotter_seqno];
854 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
858 dtc = ir->nsttcouple*ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
859 opts = &(ir->opts); /* just for ease of referencing */
862 snew(scalefac, opts->ngtc);
863 for (i = 0; i < ngtc; i++)
867 /* execute the series of trotter updates specified in the trotterpart array */
869 for (i = 0; i < NTROTTERPARTS; i++)
871 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
872 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
881 switch (trotter_seq[i])
885 boxv_trotter(ir, &(state->veta), dt, state->box, ekind, vir,
886 enerd->term[F_PDISPCORR], MassQ);
890 NHC_trotter(opts, state->nnhpres, ekind, dt, state->nhpres_xi,
891 state->nhpres_vxi, NULL, &(state->veta), MassQ, FALSE);
895 NHC_trotter(opts, opts->ngtc, ekind, dt, state->nosehoover_xi,
896 state->nosehoover_vxi, scalefac, NULL, MassQ, (ir->eI == eiVV));
897 /* need to rescale the kinetic energies and velocities here. Could
898 scale the velocities later, but we need them scaled in order to
899 produce the correct outputs, so we'll scale them here. */
901 for (t = 0; t < ngtc; t++)
903 tcstat = &ekind->tcstat[t];
904 tcstat->vscale_nhc = scalefac[t];
905 tcstat->ekinscaleh_nhc *= (scalefac[t]*scalefac[t]);
906 tcstat->ekinscalef_nhc *= (scalefac[t]*scalefac[t]);
908 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
909 /* but do we actually need the total? */
911 /* modify the velocities as well */
912 for (n = 0; n < md->homenr; n++)
914 if (md->cTC) /* does this conditional need to be here? is this always true?*/
918 for (d = 0; d < DIM; d++)
920 state->v[n][d] *= scalefac[gc];
925 for (d = 0; d < DIM; d++)
927 sumv[d] += (state->v[n][d])/md->invmass[n];
936 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
942 for (d = 0; d < DIM; d++)
944 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
946 fprintf(debug, "Conserved kappa: %15.8f %15.8f %15.8f\n", consk[0], consk[1], consk[2]);
954 extern void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit)
956 int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
957 t_grp_tcstat *tcstat;
959 real ecorr, pcorr, dvdlcorr;
960 real bmass, qmass, reft, kT, dt, ndj, nd;
961 tensor dumpres, dumvir;
963 opts = &(ir->opts); /* just for ease of referencing */
964 ngtc = ir->opts.ngtc;
965 nnhpres = state->nnhpres;
966 nh = state->nhchainlength;
972 snew(MassQ->Qinv, ngtc);
974 for (i = 0; (i < ngtc); i++)
976 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
978 MassQ->Qinv[i] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
982 MassQ->Qinv[i] = 0.0;
986 else if (EI_VV(ir->eI))
988 /* Set pressure variables */
992 if (state->vol0 == 0)
994 state->vol0 = det(state->box);
995 /* because we start by defining a fixed
996 compressibility, we need the volume at this
997 compressibility to solve the problem. */
1001 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1002 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1003 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
1004 /* An alternate mass definition, from Tuckerman et al. */
1005 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1006 for (d = 0; d < DIM; d++)
1008 for (n = 0; n < DIM; n++)
1010 MassQ->Winvm[d][n] = PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
1011 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1012 before using MTTK for anisotropic states.*/
1015 /* Allocate space for thermostat variables */
1018 snew(MassQ->Qinv, ngtc*nh);
1021 /* now, set temperature variables */
1022 for (i = 0; i < ngtc; i++)
1024 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1026 reft = max(0.0, opts->ref_t[i]);
1029 for (j = 0; j < nh; j++)
1039 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
1045 for (j = 0; j < nh; j++)
1047 MassQ->Qinv[i*nh+j] = 0.0;
1054 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
1056 int n, i, j, d, ntgrp, ngtc, nnhpres, nh, gc = 0;
1057 t_grp_tcstat *tcstat;
1059 real ecorr, pcorr, dvdlcorr;
1060 real bmass, qmass, reft, kT, dt, ndj, nd;
1061 tensor dumpres, dumvir;
1064 opts = &(ir->opts); /* just for ease of referencing */
1066 nnhpres = state->nnhpres;
1067 nh = state->nhchainlength;
1069 init_npt_masses(ir, state, MassQ, TRUE);
1071 /* first, initialize clear all the trotter calls */
1072 snew(trotter_seq, ettTSEQMAX);
1073 for (i = 0; i < ettTSEQMAX; i++)
1075 snew(trotter_seq[i], NTROTTERPARTS);
1076 for (j = 0; j < NTROTTERPARTS; j++)
1078 trotter_seq[i][j] = etrtNONE;
1080 trotter_seq[i][0] = etrtSKIPALL;
1085 /* no trotter calls, so we never use the values in the array.
1086 * We access them (so we need to define them, but ignore
1092 /* compute the kinetic energy by using the half step velocities or
1093 * the kinetic energies, depending on the order of the trotter calls */
1097 if (IR_NPT_TROTTER(ir))
1099 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1100 We start with the initial one. */
1101 /* first, a round that estimates veta. */
1102 trotter_seq[0][0] = etrtBAROV;
1104 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1106 /* The first half trotter update */
1107 trotter_seq[2][0] = etrtBAROV;
1108 trotter_seq[2][1] = etrtNHC;
1109 trotter_seq[2][2] = etrtBARONHC;
1111 /* The second half trotter update */
1112 trotter_seq[3][0] = etrtBARONHC;
1113 trotter_seq[3][1] = etrtNHC;
1114 trotter_seq[3][2] = etrtBAROV;
1116 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1119 else if (IR_NVT_TROTTER(ir))
1121 /* This is the easy version - there are only two calls, both the same.
1122 Otherwise, even easier -- no calls */
1123 trotter_seq[2][0] = etrtNHC;
1124 trotter_seq[3][0] = etrtNHC;
1126 else if (IR_NPH_TROTTER(ir))
1128 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1129 We start with the initial one. */
1130 /* first, a round that estimates veta. */
1131 trotter_seq[0][0] = etrtBAROV;
1133 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1135 /* The first half trotter update */
1136 trotter_seq[2][0] = etrtBAROV;
1137 trotter_seq[2][1] = etrtBARONHC;
1139 /* The second half trotter update */
1140 trotter_seq[3][0] = etrtBARONHC;
1141 trotter_seq[3][1] = etrtBAROV;
1143 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1146 else if (ir->eI == eiVVAK)
1148 if (IR_NPT_TROTTER(ir))
1150 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1151 We start with the initial one. */
1152 /* first, a round that estimates veta. */
1153 trotter_seq[0][0] = etrtBAROV;
1155 /* The first half trotter update, part 1 -- double update, because it commutes */
1156 trotter_seq[1][0] = etrtNHC;
1158 /* The first half trotter update, part 2 */
1159 trotter_seq[2][0] = etrtBAROV;
1160 trotter_seq[2][1] = etrtBARONHC;
1162 /* The second half trotter update, part 1 */
1163 trotter_seq[3][0] = etrtBARONHC;
1164 trotter_seq[3][1] = etrtBAROV;
1166 /* The second half trotter update */
1167 trotter_seq[4][0] = etrtNHC;
1169 else if (IR_NVT_TROTTER(ir))
1171 /* This is the easy version - there is only one call, both the same.
1172 Otherwise, even easier -- no calls */
1173 trotter_seq[1][0] = etrtNHC;
1174 trotter_seq[4][0] = etrtNHC;
1176 else if (IR_NPH_TROTTER(ir))
1178 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1179 We start with the initial one. */
1180 /* first, a round that estimates veta. */
1181 trotter_seq[0][0] = etrtBAROV;
1183 /* The first half trotter update, part 1 -- leave zero */
1184 trotter_seq[1][0] = etrtNHC;
1186 /* The first half trotter update, part 2 */
1187 trotter_seq[2][0] = etrtBAROV;
1188 trotter_seq[2][1] = etrtBARONHC;
1190 /* The second half trotter update, part 1 */
1191 trotter_seq[3][0] = etrtBARONHC;
1192 trotter_seq[3][1] = etrtBAROV;
1194 /* The second half trotter update -- blank for now */
1202 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
1205 snew(MassQ->QPinv, nnhpres*opts->nhchainlength);
1207 /* barostat temperature */
1208 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1210 reft = max(0.0, opts->ref_t[0]);
1212 for (i = 0; i < nnhpres; i++)
1214 for (j = 0; j < nh; j++)
1224 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1230 for (i = 0; i < nnhpres; i++)
1232 for (j = 0; j < nh; j++)
1234 MassQ->QPinv[i*nh+j] = 0.0;
1241 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1243 int i, j, nd, ndj, bmass, qmass, ngtcall;
1244 real ener_npt, reft, eta, kT, tau;
1247 real vol, dbaro, W, Q;
1248 int nh = state->nhchainlength;
1252 /* now we compute the contribution of the pressure to the conserved quantity*/
1254 if (ir->epc == epcMTTK)
1256 /* find the volume, and the kinetic energy of the volume */
1262 /* contribution from the pressure momenenta */
1263 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1265 /* contribution from the PV term */
1266 vol = det(state->box);
1267 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1270 case epctANISOTROPIC:
1274 case epctSURFACETENSION:
1277 case epctSEMIISOTROPIC:
1285 if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1287 /* add the energy from the barostat thermostat chain */
1288 for (i = 0; i < state->nnhpres; i++)
1291 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1292 ivxi = &state->nhpres_vxi[i*nh];
1293 ixi = &state->nhpres_xi[i*nh];
1294 iQinv = &(MassQ->QPinv[i*nh]);
1295 reft = max(ir->opts.ref_t[0], 0); /* using 'System' temperature */
1298 for (j = 0; j < nh; j++)
1302 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1303 /* contribution from the thermal variable of the NH chain */
1304 ener_npt += ixi[j]*kT;
1308 fprintf(debug, "P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f", i, j, ivxi[j], ixi[j]);
1316 for (i = 0; i < ir->opts.ngtc; i++)
1318 ixi = &state->nosehoover_xi[i*nh];
1319 ivxi = &state->nosehoover_vxi[i*nh];
1320 iQinv = &(MassQ->Qinv[i*nh]);
1322 nd = ir->opts.nrdf[i];
1323 reft = max(ir->opts.ref_t[i], 0);
1328 if (IR_NVT_TROTTER(ir))
1330 /* contribution from the thermal momenta of the NH chain */
1331 for (j = 0; j < nh; j++)
1335 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1336 /* contribution from the thermal variable of the NH chain */
1345 ener_npt += ndj*ixi[j]*kT;
1349 else /* Other non Trotter temperature NH control -- no chains yet. */
1351 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1352 ener_npt += nd*ixi[0]*kT;
1360 static real vrescale_gamdev(real ia,
1361 gmx_int64_t step, gmx_int64_t *count,
1362 gmx_int64_t seed1, gmx_int64_t seed2)
1363 /* Gamma distribution, adapted from numerical recipes */
1365 real am, e, s, v1, v2, x, y;
1376 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1378 v2 = 2.0*rnd[1] - 1.0;
1380 while (v1*v1 + v2*v2 > 1.0 ||
1381 v1*v1*GMX_REAL_MAX < 3.0*ia);
1382 /* The last check above ensures that both x (3.0 > 2.0 in s)
1383 * and the pre-factor for e do not go out of range.
1387 s = sqrt(2.0*am + 1.0);
1392 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1394 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1401 static real gaussian_count(gmx_int64_t step, gmx_int64_t *count,
1402 gmx_int64_t seed1, gmx_int64_t seed2)
1404 double rnd[2], x, y, r;
1408 gmx_rng_cycle_2uniform(step, (*count)++, seed1, seed2, rnd);
1409 x = 2.0*rnd[0] - 1.0;
1410 y = 2.0*rnd[1] - 1.0;
1413 while (r > 1.0 || r == 0.0);
1415 r = sqrt(-2.0*log(r)/r);
1420 static real vrescale_sumnoises(real nn,
1421 gmx_int64_t step, gmx_int64_t *count,
1422 gmx_int64_t seed1, gmx_int64_t seed2)
1425 * Returns the sum of nn independent gaussian noises squared
1426 * (i.e. equivalent to summing the square of the return values
1427 * of nn calls to gmx_rng_gaussian_real).
1429 const real ndeg_tol = 0.0001;
1432 if (nn < 2 + ndeg_tol)
1437 nn_int = (int)(nn + 0.5);
1439 if (nn - nn_int < -ndeg_tol || nn - nn_int > ndeg_tol)
1441 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);
1445 for (i = 0; i < nn_int; i++)
1447 gauss = gaussian_count(step, count, seed1, seed2);
1454 /* Use a gamma distribution for any real nn > 2 */
1455 r = 2.0*vrescale_gamdev(0.5*nn, step, count, seed1, seed2);
1461 static real vrescale_resamplekin(real kk, real sigma, real ndeg, real taut,
1462 gmx_int64_t step, gmx_int64_t seed)
1465 * Generates a new value for the kinetic energy,
1466 * according to Bussi et al JCP (2007), Eq. (A7)
1467 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1468 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1469 * ndeg: number of degrees of freedom of the atoms to be thermalized
1470 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1472 /* rnd_count tracks the step-local state for the cycle random
1475 gmx_int64_t rnd_count = 0;
1476 real factor, rr, ekin_new;
1480 factor = exp(-1.0/taut);
1487 rr = gaussian_count(step, &rnd_count, seed, RND_SEED_VRESCALE);
1491 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1, step, &rnd_count, seed, RND_SEED_VRESCALE) + rr*rr)/ndeg - kk) +
1492 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1497 void vrescale_tcoupl(t_inputrec *ir, gmx_int64_t step,
1498 gmx_ekindata_t *ekind, real dt,
1499 double therm_integral[])
1503 real Ek, Ek_ref1, Ek_ref, Ek_new;
1507 for (i = 0; (i < opts->ngtc); i++)
1511 Ek = trace(ekind->tcstat[i].ekinf);
1515 Ek = trace(ekind->tcstat[i].ekinh);
1518 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1520 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1521 Ek_ref = Ek_ref1*opts->nrdf[i];
1523 Ek_new = vrescale_resamplekin(Ek, Ek_ref, opts->nrdf[i],
1527 /* Analytically Ek_new>=0, but we check for rounding errors */
1530 ekind->tcstat[i].lambda = 0.0;
1534 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1537 therm_integral[i] -= Ek_new - Ek;
1541 fprintf(debug, "TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1542 i, Ek_ref, Ek, Ek_new, ekind->tcstat[i].lambda);
1547 ekind->tcstat[i].lambda = 1.0;
1552 real vrescale_energy(t_grpopts *opts, double therm_integral[])
1558 for (i = 0; i < opts->ngtc; i++)
1560 ener += therm_integral[i];
1566 void rescale_velocities(gmx_ekindata_t *ekind, t_mdatoms *mdatoms,
1567 int start, int end, rvec v[])
1570 t_grp_tcstat *tcstat;
1571 unsigned short *cACC, *cTC;
1576 tcstat = ekind->tcstat;
1581 gstat = ekind->grpstat;
1582 cACC = mdatoms->cACC;
1586 for (n = start; n < end; n++)
1596 /* Only scale the velocity component relative to the COM velocity */
1597 rvec_sub(v[n], gstat[ga].u, vrel);
1598 lg = tcstat[gt].lambda;
1599 for (d = 0; d < DIM; d++)
1601 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1608 for (n = start; n < end; n++)
1614 lg = tcstat[gt].lambda;
1615 for (d = 0; d < DIM; d++)
1624 /* set target temperatures if we are annealing */
1625 void update_annealing_target_temp(t_grpopts *opts, real t)
1627 int i, j, n, npoints;
1628 real pert, thist = 0, x;
1630 for (i = 0; i < opts->ngtc; i++)
1632 npoints = opts->anneal_npoints[i];
1633 switch (opts->annealing[i])
1638 /* calculate time modulo the period */
1639 pert = opts->anneal_time[i][npoints-1];
1641 thist = t - n*pert; /* modulo time */
1642 /* Make sure rounding didn't get us outside the interval */
1643 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1652 gmx_fatal(FARGS, "Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)", i, opts->ngtc, npoints);
1654 /* We are doing annealing for this group if we got here,
1655 * and we have the (relative) time as thist.
1656 * calculate target temp */
1658 while ((j < npoints-1) && (thist > (opts->anneal_time[i][j+1])))
1664 /* Found our position between points j and j+1.
1665 * Interpolate: x is the amount from j+1, (1-x) from point j
1666 * First treat possible jumps in temperature as a special case.
1668 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1670 opts->ref_t[i] = opts->anneal_temp[i][j+1];
1674 x = ((thist-opts->anneal_time[i][j])/
1675 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1676 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1681 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];