1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * GROwing Monsters And Cloning Shrimps
47 #include "gmx_fatal.h"
50 #include "gmx_random.h"
54 #define NTROTTERPARTS 3
56 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
58 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
59 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
61 #define MAX_SUZUKI_YOSHIDA_NUM 5
62 #define SUZUKI_YOSHIDA_NUM 5
64 static const double sy_const_1[] = { 1. };
65 static const double sy_const_3[] = { 0.828981543588751,-0.657963087177502,0.828981543588751 };
66 static const double sy_const_5[] = { 0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065 };
68 static const double* sy_const[] = {
78 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
82 {0.828981543588751,-0.657963087177502,0.828981543588751},
84 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
87 /* these integration routines are only referenced inside this file */
88 static void NHC_trotter(t_grpopts *opts,int nvar, gmx_ekindata_t *ekind,real dtfull,
89 double xi[],double vxi[], double scalefac[], real *veta, t_extmass *MassQ, gmx_bool bEkinAveVel)
92 /* general routine for both barostat and thermostat nose hoover chains */
95 double Ekin,Efac,reft,kT,nd;
102 int mstepsi, mstepsj;
103 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
104 int nh = opts->nhchainlength;
107 mstepsi = mstepsj = ns;
109 /* if scalefac is NULL, we are doing the NHC of the barostat */
112 if (scalefac == NULL) {
116 for (i=0; i<nvar; i++)
119 /* make it easier to iterate by selecting
120 out the sub-array that corresponds to this T group */
125 iQinv = &(MassQ->QPinv[i*nh]);
126 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
127 reft = max(0.0,opts->ref_t[0]);
128 Ekin = sqr(*veta)/MassQ->Winv;
130 iQinv = &(MassQ->Qinv[i*nh]);
131 tcstat = &ekind->tcstat[i];
133 reft = max(0.0,opts->ref_t[i]);
136 Ekin = 2*trace(tcstat->ekinf)*tcstat->ekinscalef_nhc;
138 Ekin = 2*trace(tcstat->ekinh)*tcstat->ekinscaleh_nhc;
143 for(mi=0;mi<mstepsi;mi++)
145 for(mj=0;mj<mstepsj;mj++)
147 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
148 dt = sy_const[ns][mj] * dtfull / mstepsi;
150 /* compute the thermal forces */
151 GQ[0] = iQinv[0]*(Ekin - nd*kT);
155 if (iQinv[j+1] > 0) {
156 /* we actually don't need to update here if we save the
157 state of the GQ, but it's easier to just recompute*/
158 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
164 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
167 Efac = exp(-0.125*dt*ivxi[j]);
168 ivxi[j-1] = Efac*(ivxi[j-1]*Efac + 0.25*dt*GQ[j-1]);
171 Efac = exp(-0.5*dt*ivxi[0]);
179 /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
180 able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
181 think about this a bit more . . . */
183 GQ[0] = iQinv[0]*(Ekin - nd*kT);
185 /* update thermostat positions */
188 ixi[j] += 0.5*dt*ivxi[j];
193 Efac = exp(-0.125*dt*ivxi[j+1]);
194 ivxi[j] = Efac*(ivxi[j]*Efac + 0.25*dt*GQ[j]);
195 if (iQinv[j+1] > 0) {
196 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
201 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
208 static void boxv_trotter(t_inputrec *ir, real *veta, real dt, tensor box,
209 gmx_ekindata_t *ekind, tensor vir, real pcorr, real ecorr, t_extmass *MassQ)
216 tensor Winvm,ekinmod,localpres;
218 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
219 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
222 if (ir->epct==epctSEMIISOTROPIC)
231 /* eta is in pure units. veta is in units of ps^-1. GW is in
232 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
233 taken to use only RATIOS of eta in updating the volume. */
235 /* we take the partial pressure tensors, modify the
236 kinetic energy tensor, and recovert to pressure */
238 if (ir->opts.nrdf[0]==0)
240 gmx_fatal(FARGS,"Barostat is coupled to a T-group with no degrees of freedom\n");
242 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
243 alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
244 alpha *= ekind->tcstat[0].ekinscalef_nhc;
245 msmul(ekind->ekin,alpha,ekinmod);
246 /* for now, we use Elr = 0, because if you want to get it right, you
247 really should be using PME. Maybe print a warning? */
249 pscal = calc_pres(ir->ePBC,nwall,box,ekinmod,vir,localpres);
252 GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
258 * This file implements temperature and pressure coupling algorithms:
259 * For now only the Weak coupling and the modified weak coupling.
261 * Furthermore computation of pressure and temperature is done here
265 real calc_pres(int ePBC,int nwall,matrix box,tensor ekin,tensor vir,
271 if (ePBC==epbcNONE || (ePBC==epbcXY && nwall!=2))
274 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
275 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
279 fac=PRESFAC*2.0/det(box);
280 for(n=0; (n<DIM); n++)
281 for(m=0; (m<DIM); m++)
282 pres[n][m] = (ekin[n][m] - vir[n][m])*fac;
285 pr_rvecs(debug,0,"PC: pres",pres,DIM);
286 pr_rvecs(debug,0,"PC: ekin",ekin,DIM);
287 pr_rvecs(debug,0,"PC: vir ",vir, DIM);
288 pr_rvecs(debug,0,"PC: box ",box, DIM);
291 return trace(pres)/DIM;
294 real calc_temp(real ekin,real nrdf)
297 return (2.0*ekin)/(nrdf*BOLTZ);
302 void parrinellorahman_pcoupl(FILE *fplog,gmx_large_int_t step,
303 t_inputrec *ir,real dt,tensor pres,
304 tensor box,tensor box_rel,tensor boxv,
305 tensor M,matrix mu,gmx_bool bFirstStep)
307 /* This doesn't do any coordinate updating. It just
308 * integrates the box vector equations from the calculated
309 * acceleration due to pressure difference. We also compute
310 * the tensor M which is used in update to couple the particle
311 * coordinates to the box vectors.
313 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
316 * M_nk = (h') * (h' * h + h' h) * h
318 * with the dots denoting time derivatives and h is the transformation from
319 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
320 * This also goes for the pressure and M tensors - they are transposed relative
321 * to ours. Our equation thus becomes:
324 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
326 * where b is the gromacs box matrix.
327 * Our box accelerations are given by
329 * b = vol/W inv(box') * (P-ref_P) (=h')
334 real vol=box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
335 real atot,arel,change,maxchange,xy_pressure;
336 tensor invbox,pdiff,t1,t2;
340 m_inv_ur0(box,invbox);
343 /* Note that PRESFAC does not occur here.
344 * The pressure and compressibility always occur as a product,
345 * therefore the pressure unit drops out.
347 maxl=max(box[XX][XX],box[YY][YY]);
348 maxl=max(maxl,box[ZZ][ZZ]);
352 (4*M_PI*M_PI*ir->compress[d][n])/(3*ir->tau_p*ir->tau_p*maxl);
354 m_sub(pres,ir->ref_p,pdiff);
356 if(ir->epct==epctSURFACETENSION) {
357 /* Unlike Berendsen coupling it might not be trivial to include a z
358 * pressure correction here? On the other hand we don't scale the
359 * box momentarily, but change accelerations, so it might not be crucial.
361 xy_pressure=0.5*(pres[XX][XX]+pres[YY][YY]);
363 pdiff[d][d]=(xy_pressure-(pres[ZZ][ZZ]-ir->ref_p[d][d]/box[d][d]));
366 tmmul(invbox,pdiff,t1);
367 /* Move the off-diagonal elements of the 'force' to one side to ensure
368 * that we obey the box constraints.
372 t1[d][n] += t1[n][d];
378 case epctANISOTROPIC:
381 t1[d][n] *= winv[d][n]*vol;
384 /* calculate total volume acceleration */
385 atot=box[XX][XX]*box[YY][YY]*t1[ZZ][ZZ]+
386 box[XX][XX]*t1[YY][YY]*box[ZZ][ZZ]+
387 t1[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
389 /* set all RELATIVE box accelerations equal, and maintain total V
393 t1[d][n] = winv[0][0]*vol*arel*box[d][n];
395 case epctSEMIISOTROPIC:
396 case epctSURFACETENSION:
397 /* Note the correction to pdiff above for surftens. coupling */
399 /* calculate total XY volume acceleration */
400 atot=box[XX][XX]*t1[YY][YY]+t1[XX][XX]*box[YY][YY];
401 arel=atot/(2*box[XX][XX]*box[YY][YY]);
402 /* set RELATIVE XY box accelerations equal, and maintain total V
403 * change speed. Dont change the third box vector accelerations */
406 t1[d][n] = winv[d][n]*vol*arel*box[d][n];
408 t1[ZZ][n] *= winv[d][n]*vol;
411 gmx_fatal(FARGS,"Parrinello-Rahman pressure coupling type %s "
412 "not supported yet\n",EPCOUPLTYPETYPE(ir->epct));
419 boxv[d][n] += dt*t1[d][n];
421 /* We do NOT update the box vectors themselves here, since
422 * we need them for shifting later. It is instead done last
423 * in the update() routine.
426 /* Calculate the change relative to diagonal elements-
427 since it's perfectly ok for the off-diagonal ones to
428 be zero it doesn't make sense to check the change relative
432 change=fabs(dt*boxv[d][n]/box[d][d]);
434 if (change>maxchange)
438 if (maxchange > 0.01 && fplog) {
441 "\nStep %s Warning: Pressure scaling more than 1%%. "
442 "This may mean your system\n is not yet equilibrated. "
443 "Use of Parrinello-Rahman pressure coupling during\n"
444 "equilibration can lead to simulation instability, "
445 "and is discouraged.\n",
446 gmx_step_str(step,buf));
450 preserve_box_shape(ir,box_rel,boxv);
452 mtmul(boxv,box,t1); /* t1=boxv * b' */
456 /* Determine the scaling matrix mu for the coordinates */
459 t1[d][n] = box[d][n] + dt*boxv[d][n];
460 preserve_box_shape(ir,box_rel,t1);
461 /* t1 is the box at t+dt, determine mu as the relative change */
462 mmul_ur0(invbox,t1,mu);
465 void berendsen_pcoupl(FILE *fplog,gmx_large_int_t step,
466 t_inputrec *ir,real dt, tensor pres,matrix box,
470 real scalar_pressure, xy_pressure, p_corr_z;
471 char *ptr,buf[STRLEN];
474 * Calculate the scaling matrix mu
478 for(d=0; d<DIM; d++) {
479 scalar_pressure += pres[d][d]/DIM;
481 xy_pressure += pres[d][d]/(DIM-1);
483 /* Pressure is now in bar, everywhere. */
484 #define factor(d,m) (ir->compress[d][m]*dt/ir->tau_p)
486 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
487 * necessary for triclinic scaling
494 mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
497 case epctSEMIISOTROPIC:
499 mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
501 1.0 - factor(ZZ,ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
503 case epctANISOTROPIC:
506 mu[d][n] = (d==n ? 1.0 : 0.0)
507 -factor(d,n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
509 case epctSURFACETENSION:
510 /* ir->ref_p[0/1] is the reference surface-tension times *
511 * the number of surfaces */
512 if (ir->compress[ZZ][ZZ])
513 p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
515 /* when the compressibity is zero, set the pressure correction *
516 * in the z-direction to zero to get the correct surface tension */
518 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
519 for(d=0; d<DIM-1; d++)
520 mu[d][d] = 1.0 + factor(d,d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
521 - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
524 gmx_fatal(FARGS,"Berendsen pressure coupling type %s not supported yet\n",
525 EPCOUPLTYPETYPE(ir->epct));
528 /* To fullfill the orientation restrictions on triclinic boxes
529 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
530 * the other elements of mu to first order.
532 mu[YY][XX] += mu[XX][YY];
533 mu[ZZ][XX] += mu[XX][ZZ];
534 mu[ZZ][YY] += mu[YY][ZZ];
540 pr_rvecs(debug,0,"PC: pres ",pres,3);
541 pr_rvecs(debug,0,"PC: mu ",mu,3);
544 if (mu[XX][XX]<0.99 || mu[XX][XX]>1.01 ||
545 mu[YY][YY]<0.99 || mu[YY][YY]>1.01 ||
546 mu[ZZ][ZZ]<0.99 || mu[ZZ][ZZ]>1.01) {
548 sprintf(buf,"\nStep %s Warning: pressure scaling more than 1%%, "
550 gmx_step_str(step,buf2),mu[XX][XX],mu[YY][YY],mu[ZZ][ZZ]);
552 fprintf(fplog,"%s",buf);
553 fprintf(stderr,"%s",buf);
557 void berendsen_pscale(t_inputrec *ir,matrix mu,
558 matrix box,matrix box_rel,
559 int start,int nr_atoms,
560 rvec x[],unsigned short cFREEZE[],
563 ivec *nFreeze=ir->opts.nFreeze;
566 /* Scale the positions */
567 for (n=start; n<start+nr_atoms; n++) {
572 x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
574 x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
576 x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
578 /* compute final boxlengths */
579 for (d=0; d<DIM; d++) {
580 box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
581 box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
582 box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
585 preserve_box_shape(ir,box_rel,box);
587 /* (un)shifting should NOT be done after this,
588 * since the box vectors might have changed
590 inc_nrnb(nrnb,eNR_PCOUPL,nr_atoms);
593 void berendsen_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt)
601 for(i=0; (i<opts->ngtc); i++)
605 T = ekind->tcstat[i].T;
609 T = ekind->tcstat[i].Th;
612 if ((opts->tau_t[i] > 0) && (T > 0.0)) {
613 reft = max(0.0,opts->ref_t[i]);
614 lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
615 ekind->tcstat[i].lambda = max(min(lll,1.25),0.8);
618 ekind->tcstat[i].lambda = 1.0;
623 fprintf(debug,"TC: group %d: T: %g, Lambda: %g\n",
624 i,T,ekind->tcstat[i].lambda);
629 static int poisson_variate(real lambda,gmx_rng_t rng) {
640 p *= gmx_rng_uniform_real(rng);
646 void andersen_tcoupl(t_inputrec *ir,t_mdatoms *md,t_state *state, gmx_rng_t rng, real rate, t_idef *idef, int nblocks, int *sblock,gmx_bool *randatom, int *randatom_list, gmx_bool *randomize, real *boltzfac)
649 int i,j,k,d,len,n,ngtc,gc=0;
650 int nshake, nsettle, nrandom, nrand_group;
651 real boltz,scal,reft,prand;
654 /* convenience variables */
658 /* idef is only passed in if it's chance-per-particle andersen, so
659 it essentially serves as a boolean to determine which type of
660 andersen is being used */
663 /* randomly atoms to randomize. However, all constraint
664 groups have to have either all of the atoms or none of the
668 1. Select whether or not to randomize each atom to get the correct probability.
669 2. Cycle through the constraint groups.
670 2a. for each constraint group, determine the fraction f of that constraint group that are
671 chosen to be randomized.
672 2b. all atoms in the constraint group are randomized with probability f.
676 if ((rate < 0.05) && (md->homenr > 50))
678 /* if the rate is relatively high, use a standard method, if low rate,
680 /* poisson distributions approxmation, more efficient for
681 * low rates, fewer random numbers required */
682 nrandom = poisson_variate(md->homenr*rate,rng); /* how many do we randomize? Use Poisson. */
683 /* now we know how many, choose them randomly. No worries about repeats, at this rate, it's negligible.
684 worst thing that happens, it lowers the true rate an negligible amount */
685 for (i=0;i<nrandom;i++)
687 randatom[(int)(gmx_rng_uniform_real(rng)*md->homenr)] = TRUE;
692 for (i=0;i<md->homenr;i++)
694 if (gmx_rng_uniform_real(rng)<rate)
702 /* instead of looping over the constraint groups, if we had a
703 list of which atoms were in which constraint groups, we
704 could then loop over only the groups that are randomized
705 now. But that is not available now. Create later after
706 determining whether there actually is any slowing. */
708 /* first, loop through the settles to make sure all groups either entirely randomized, or not randomized. */
710 nsettle = idef->il[F_SETTLE].nr/2;
711 for (i=0;i<nsettle;i++)
713 iatoms = idef->il[F_SETTLE].iatoms;
715 for (k=0;k<3;k++) /* settles are always 3 atoms, hardcoded */
717 if (randatom[iatoms[2*i+1]+k])
719 nrand_group++; /* count the number of atoms to be shaken in the settles group */
720 randatom[iatoms[2*i+1]+k] = FALSE;
726 prand = (nrand_group)/3.0; /* use this fraction to compute the probability the
727 whole group is randomized */
728 if (gmx_rng_uniform_real(rng)<prand)
732 randatom[iatoms[2*i+1]+k] = TRUE; /* mark them all to be randomized */
739 /* now loop through the shake groups */
741 for (i=0;i<nshake;i++)
743 iatoms = &(idef->il[F_CONSTR].iatoms[sblock[i]]);
744 len = sblock[i+1]-sblock[i];
749 { /* only 2/3 of the sblock items are atoms, the others are labels */
750 if (randatom[iatoms[k]])
753 randatom[iatoms[k]] = FALSE; /* need to mark it false here in case the atom is in more than
754 one group in the shake block */
761 prand = (nrand_group)/(1.0*(2*len/3));
762 if (gmx_rng_uniform_real(rng)<prand)
767 { /* only 2/3 of the sblock items are atoms, the others are labels */
768 randatom[iatoms[k]] = TRUE; /* randomize all of them */
778 for (i=0;i<md->homenr;i++) /* now loop over the list of atoms */
782 randatom_list[n] = i;
786 nrandom = n; /* there are some values of nrandom for which
787 this algorithm won't work; for example all
788 water molecules and nrandom =/= 3. Better to
789 recount and use this number (which we
790 calculate anyway: it will not affect
791 the average number of atoms accepted.
797 /* if it's andersen-massive, then randomize all the atoms */
798 nrandom = md->homenr;
799 for (i=0;i<nrandom;i++)
801 randatom_list[i] = i;
805 /* randomize the velocities of the selected particles */
807 for (i=0;i<nrandom;i++) /* now loop over the list of atoms */
809 n = randatom_list[i];
812 gc = md->cTC[n]; /* assign the atom to a temperature group if there are more than one */
816 scal = sqrt(boltzfac[gc]*md->invmass[n]);
819 state->v[n][d] = scal*gmx_rng_gaussian_table(rng);
822 randatom[n] = FALSE; /* unmark this atom for randomization */
827 void nosehoover_tcoupl(t_grpopts *opts,gmx_ekindata_t *ekind,real dt,
828 double xi[],double vxi[], t_extmass *MassQ)
833 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
835 for(i=0; (i<opts->ngtc); i++)
837 reft = max(0.0,opts->ref_t[i]);
839 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
840 xi[i] += dt*(oldvxi + vxi[i])*0.5;
844 t_state *init_bufstate(const t_state *template_state)
847 int nc = template_state->nhchainlength;
849 snew(state->nosehoover_xi,nc*template_state->ngtc);
850 snew(state->nosehoover_vxi,nc*template_state->ngtc);
851 snew(state->therm_integral,template_state->ngtc);
852 snew(state->nhpres_xi,nc*template_state->nnhpres);
853 snew(state->nhpres_vxi,nc*template_state->nnhpres);
858 void destroy_bufstate(t_state *state)
862 sfree(state->nosehoover_xi);
863 sfree(state->nosehoover_vxi);
864 sfree(state->therm_integral);
865 sfree(state->nhpres_xi);
866 sfree(state->nhpres_vxi);
870 void trotter_update(t_inputrec *ir,gmx_large_int_t step, gmx_ekindata_t *ekind,
871 gmx_enerdata_t *enerd, t_state *state,
872 tensor vir, t_mdatoms *md,
873 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
876 int n,i,j,d,ntgrp,ngtc,gc=0;
877 t_grp_tcstat *tcstat;
879 gmx_large_int_t step_eff;
880 real ecorr,pcorr,dvdlcorr;
881 real bmass,qmass,reft,kT,dt,nd;
882 tensor dumpres,dumvir;
883 double *scalefac,dtc;
885 rvec sumv={0,0,0},consk;
888 if (trotter_seqno <= ettTSEQ2)
890 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
891 is actually the last half step from the previous step. Thus the first half step
892 actually corresponds to the n-1 step*/
898 bCouple = (ir->nsttcouple == 1 ||
899 do_per_step(step_eff+ir->nsttcouple,ir->nsttcouple));
901 trotter_seq = trotter_seqlist[trotter_seqno];
903 /* signal we are returning if nothing is going to be done in this routine */
904 if ((trotter_seq[0] == etrtSKIPALL) || !(bCouple))
909 dtc = ir->nsttcouple*ir->delta_t;
910 opts = &(ir->opts); /* just for ease of referencing */
913 snew(scalefac,opts->ngtc);
918 /* execute the series of trotter updates specified in the trotterpart array */
920 for (i=0;i<NTROTTERPARTS;i++){
921 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
922 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
931 switch (trotter_seq[i])
935 boxv_trotter(ir,&(state->veta),dt,state->box,ekind,vir,
936 enerd->term[F_PDISPCORR],enerd->term[F_DISPCORR],MassQ);
940 NHC_trotter(opts,state->nnhpres,ekind,dt,state->nhpres_xi,
941 state->nhpres_vxi,NULL,&(state->veta),MassQ,FALSE);
945 NHC_trotter(opts,opts->ngtc,ekind,dt,state->nosehoover_xi,
946 state->nosehoover_vxi,scalefac,NULL,MassQ,(ir->eI==eiVV));
947 /* need to rescale the kinetic energies and velocities here. Could
948 scale the velocities later, but we need them scaled in order to
949 produce the correct outputs, so we'll scale them here. */
951 for (i=0; i<ngtc;i++)
953 tcstat = &ekind->tcstat[i];
954 tcstat->vscale_nhc = scalefac[i];
955 tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
956 tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
958 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
959 /* but do we actually need the total? */
961 /* modify the velocities as well */
962 for (n=md->start;n<md->start+md->homenr;n++)
964 if (md->cTC) /* does this conditional need to be here? is this always true?*/
970 state->v[n][d] *= scalefac[gc];
977 sumv[d] += (state->v[n][d])/md->invmass[n];
986 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
994 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
996 fprintf(debug,"Conserved kappa: %15.8f %15.8f %15.8f\n",consk[0],consk[1],consk[2]);
1004 extern void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit)
1006 int n,i,j,d,ntgrp,ngtc,nnhpres,nh,gc=0;
1007 t_grp_tcstat *tcstat;
1009 real ecorr,pcorr,dvdlcorr;
1010 real bmass,qmass,reft,kT,dt,ndj,nd;
1011 tensor dumpres,dumvir;
1013 opts = &(ir->opts); /* just for ease of referencing */
1014 ngtc = ir->opts.ngtc;
1015 nnhpres = state->nnhpres;
1016 nh = state->nhchainlength;
1018 if (ir->eI == eiMD) {
1021 snew(MassQ->Qinv,ngtc);
1023 for(i=0; (i<ngtc); i++)
1025 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1027 MassQ->Qinv[i]=1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
1035 else if (EI_VV(ir->eI))
1037 /* Set pressure variables */
1041 if (state->vol0 == 0)
1043 state->vol0 = det(state->box);
1044 /* because we start by defining a fixed
1045 compressibility, we need the volume at this
1046 compressibility to solve the problem. */
1050 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1051 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1052 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
1053 /* An alternate mass definition, from Tuckerman et al. */
1054 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1059 MassQ->Winvm[d][n]= PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
1060 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1061 before using MTTK for anisotropic states.*/
1064 /* Allocate space for thermostat variables */
1067 snew(MassQ->Qinv,ngtc*nh);
1070 /* now, set temperature variables */
1071 for (i=0; i<ngtc; i++)
1073 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1075 reft = max(0.0,opts->ref_t[i]);
1088 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
1096 MassQ->Qinv[i*nh+j] = 0.0;
1103 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
1105 int n,i,j,d,ntgrp,ngtc,nnhpres,nh,gc=0;
1106 t_grp_tcstat *tcstat;
1108 real ecorr,pcorr,dvdlcorr;
1109 real bmass,qmass,reft,kT,dt,ndj,nd;
1110 tensor dumpres,dumvir;
1113 opts = &(ir->opts); /* just for ease of referencing */
1115 nnhpres = state->nnhpres;
1116 nh = state->nhchainlength;
1118 init_npt_masses(ir, state, MassQ, TRUE);
1120 /* first, initialize clear all the trotter calls */
1121 snew(trotter_seq,ettTSEQMAX);
1122 for (i=0;i<ettTSEQMAX;i++)
1124 snew(trotter_seq[i],NTROTTERPARTS);
1125 for (j=0;j<NTROTTERPARTS;j++) {
1126 trotter_seq[i][j] = etrtNONE;
1128 trotter_seq[i][0] = etrtSKIPALL;
1133 /* no trotter calls, so we never use the values in the array.
1134 * We access them (so we need to define them, but ignore
1140 /* compute the kinetic energy by using the half step velocities or
1141 * the kinetic energies, depending on the order of the trotter calls */
1145 if (IR_NPT_TROTTER(ir))
1147 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1148 We start with the initial one. */
1149 /* first, a round that estimates veta. */
1150 trotter_seq[0][0] = etrtBAROV;
1152 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1154 /* The first half trotter update */
1155 trotter_seq[2][0] = etrtBAROV;
1156 trotter_seq[2][1] = etrtNHC;
1157 trotter_seq[2][2] = etrtBARONHC;
1159 /* The second half trotter update */
1160 trotter_seq[3][0] = etrtBARONHC;
1161 trotter_seq[3][1] = etrtNHC;
1162 trotter_seq[3][2] = etrtBAROV;
1164 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1167 else if (IR_NVT_TROTTER(ir))
1169 /* This is the easy version - there are only two calls, both the same.
1170 Otherwise, even easier -- no calls */
1171 trotter_seq[2][0] = etrtNHC;
1172 trotter_seq[3][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 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1183 /* The first half trotter update */
1184 trotter_seq[2][0] = etrtBAROV;
1185 trotter_seq[2][1] = etrtBARONHC;
1187 /* The second half trotter update */
1188 trotter_seq[3][0] = etrtBARONHC;
1189 trotter_seq[3][1] = etrtBAROV;
1191 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1194 else if (ir->eI==eiVVAK)
1196 if (IR_NPT_TROTTER(ir))
1198 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1199 We start with the initial one. */
1200 /* first, a round that estimates veta. */
1201 trotter_seq[0][0] = etrtBAROV;
1203 /* The first half trotter update, part 1 -- double update, because it commutes */
1204 trotter_seq[1][0] = etrtNHC;
1206 /* The first half trotter update, part 2 */
1207 trotter_seq[2][0] = etrtBAROV;
1208 trotter_seq[2][1] = etrtBARONHC;
1210 /* The second half trotter update, part 1 */
1211 trotter_seq[3][0] = etrtBARONHC;
1212 trotter_seq[3][1] = etrtBAROV;
1214 /* The second half trotter update */
1215 trotter_seq[4][0] = etrtNHC;
1217 else if (IR_NVT_TROTTER(ir))
1219 /* This is the easy version - there is only one call, both the same.
1220 Otherwise, even easier -- no calls */
1221 trotter_seq[1][0] = etrtNHC;
1222 trotter_seq[4][0] = etrtNHC;
1224 else if (IR_NPH_TROTTER(ir))
1226 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1227 We start with the initial one. */
1228 /* first, a round that estimates veta. */
1229 trotter_seq[0][0] = etrtBAROV;
1231 /* The first half trotter update, part 1 -- leave zero */
1232 trotter_seq[1][0] = etrtNHC;
1234 /* The first half trotter update, part 2 */
1235 trotter_seq[2][0] = etrtBAROV;
1236 trotter_seq[2][1] = etrtBARONHC;
1238 /* The second half trotter update, part 1 */
1239 trotter_seq[3][0] = etrtBARONHC;
1240 trotter_seq[3][1] = etrtBAROV;
1242 /* The second half trotter update -- blank for now */
1250 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
1253 snew(MassQ->QPinv,nnhpres*opts->nhchainlength);
1255 /* barostat temperature */
1256 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1258 reft = max(0.0,opts->ref_t[0]);
1260 for (i=0;i<nnhpres;i++) {
1270 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1276 for (i=0;i<nnhpres;i++) {
1279 MassQ->QPinv[i*nh+j]=0.0;
1286 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1288 int i,j,nd,ndj,bmass,qmass,ngtcall;
1289 real ener_npt,reft,eta,kT,tau;
1293 int nh = state->nhchainlength;
1297 /* now we compute the contribution of the pressure to the conserved quantity*/
1299 if (ir->epc==epcMTTK)
1301 /* find the volume, and the kinetic energy of the volume */
1306 /* contribution from the pressure momenenta */
1307 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1309 /* contribution from the PV term */
1310 vol = det(state->box);
1311 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1314 case epctANISOTROPIC:
1318 case epctSURFACETENSION:
1321 case epctSEMIISOTROPIC:
1329 if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1331 /* add the energy from the barostat thermostat chain */
1332 for (i=0;i<state->nnhpres;i++) {
1334 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1335 ivxi = &state->nhpres_vxi[i*nh];
1336 ixi = &state->nhpres_xi[i*nh];
1337 iQinv = &(MassQ->QPinv[i*nh]);
1338 reft = max(ir->opts.ref_t[0],0); /* using 'System' temperature */
1345 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1346 /* contribution from the thermal variable of the NH chain */
1347 ener_npt += ixi[j]*kT;
1351 fprintf(debug,"P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f",i,j,ivxi[j],ixi[j]);
1359 for(i=0; i<ir->opts.ngtc; i++)
1361 ixi = &state->nosehoover_xi[i*nh];
1362 ivxi = &state->nosehoover_vxi[i*nh];
1363 iQinv = &(MassQ->Qinv[i*nh]);
1365 nd = ir->opts.nrdf[i];
1366 reft = max(ir->opts.ref_t[i],0);
1371 if (IR_NVT_TROTTER(ir))
1373 /* contribution from the thermal momenta of the NH chain */
1377 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1378 /* contribution from the thermal variable of the NH chain */
1386 ener_npt += ndj*ixi[j]*kT;
1390 else /* Other non Trotter temperature NH control -- no chains yet. */
1392 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1393 ener_npt += nd*ixi[0]*kT;
1401 static real vrescale_gamdev(int ia, gmx_rng_t rng)
1402 /* Gamma distribution, adapted from numerical recipes */
1405 real am,e,s,v1,v2,x,y;
1412 for(j=1; j<=ia; j++)
1414 x *= gmx_rng_uniform_real(rng);
1428 v1 = gmx_rng_uniform_real(rng);
1429 v2 = 2.0*gmx_rng_uniform_real(rng)-1.0;
1431 while (v1*v1 + v2*v2 > 1.0 ||
1432 v1*v1*GMX_REAL_MAX < 3.0*ia);
1433 /* The last check above ensures that both x (3.0 > 2.0 in s)
1434 * and the pre-factor for e do not go out of range.
1438 s = sqrt(2.0*am + 1.0);
1442 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1444 while (gmx_rng_uniform_real(rng) > e);
1450 static real vrescale_sumnoises(int nn,gmx_rng_t rng)
1453 * Returns the sum of n independent gaussian noises squared
1454 * (i.e. equivalent to summing the square of the return values
1455 * of nn calls to gmx_rng_gaussian_real).xs
1461 } else if (nn == 1) {
1462 rr = gmx_rng_gaussian_real(rng);
1464 } else if (nn % 2 == 0) {
1465 return 2.0*vrescale_gamdev(nn/2,rng);
1467 rr = gmx_rng_gaussian_real(rng);
1468 return 2.0*vrescale_gamdev((nn-1)/2,rng) + rr*rr;
1472 static real vrescale_resamplekin(real kk,real sigma, int ndeg, real taut,
1476 * Generates a new value for the kinetic energy,
1477 * according to Bussi et al JCP (2007), Eq. (A7)
1478 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1479 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1480 * ndeg: number of degrees of freedom of the atoms to be thermalized
1481 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1486 factor = exp(-1.0/taut);
1490 rr = gmx_rng_gaussian_real(rng);
1493 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1,rng) + rr*rr)/ndeg - kk) +
1494 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1497 void vrescale_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt,
1498 double therm_integral[],gmx_rng_t rng)
1502 real Ek,Ek_ref1,Ek_ref,Ek_new;
1506 for(i=0; (i<opts->ngtc); i++)
1510 Ek = trace(ekind->tcstat[i].ekinf);
1514 Ek = trace(ekind->tcstat[i].ekinh);
1517 if (opts->tau_t[i] > 0 && opts->nrdf[i] > 0 && Ek > 0)
1519 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1520 Ek_ref = Ek_ref1*opts->nrdf[i];
1522 Ek_new = vrescale_resamplekin(Ek,Ek_ref,opts->nrdf[i],
1523 opts->tau_t[i]/dt,rng);
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++) {
1557 ener += therm_integral[i];
1563 void rescale_velocities(gmx_ekindata_t *ekind,t_mdatoms *mdatoms,
1564 int start,int end,rvec v[])
1567 t_grp_tcstat *tcstat;
1568 unsigned short *cACC,*cTC;
1573 tcstat = ekind->tcstat;
1578 gstat = ekind->grpstat;
1579 cACC = mdatoms->cACC;
1583 for(n=start; n<end; n++)
1593 /* Only scale the velocity component relative to the COM velocity */
1594 rvec_sub(v[n],gstat[ga].u,vrel);
1595 lg = tcstat[gt].lambda;
1596 for(d=0; d<DIM; d++)
1598 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1605 for(n=start; n<end; n++)
1611 lg = tcstat[gt].lambda;
1612 for(d=0; d<DIM; d++)
1621 /* set target temperatures if we are annealing */
1622 void update_annealing_target_temp(t_grpopts *opts,real t)
1625 real pert,thist=0,x;
1627 for(i=0;i<opts->ngtc;i++) {
1628 npoints = opts->anneal_npoints[i];
1629 switch (opts->annealing[i]) {
1633 /* calculate time modulo the period */
1634 pert = opts->anneal_time[i][npoints-1];
1636 thist = t - n*pert; /* modulo time */
1637 /* Make sure rounding didn't get us outside the interval */
1638 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1645 gmx_fatal(FARGS,"Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",i,opts->ngtc,npoints);
1647 /* We are doing annealing for this group if we got here,
1648 * and we have the (relative) time as thist.
1649 * calculate target temp */
1651 while ((j < npoints-1) && (thist>(opts->anneal_time[i][j+1])))
1653 if (j < npoints-1) {
1654 /* Found our position between points j and j+1.
1655 * Interpolate: x is the amount from j+1, (1-x) from point j
1656 * First treat possible jumps in temperature as a special case.
1658 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1659 opts->ref_t[i]=opts->anneal_temp[i][j+1];
1661 x = ((thist-opts->anneal_time[i][j])/
1662 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1663 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1667 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];