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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
50 #include "gmx_fatal.h"
53 #include "gmx_random.h"
57 #define NTROTTERPARTS 3
59 /* these integration routines are only referenced inside this file */
60 static void NHC_trotter(t_grpopts *opts,int nvar, gmx_ekindata_t *ekind,real dtfull,
61 double xi[],double vxi[], double scalefac[], real *veta, t_extmass *MassQ, gmx_bool bEkinAveVel)
64 /* general routine for both barostat and thermostat nose hoover chains */
67 double Ekin,Efac,reft,kT,nd;
75 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
76 int nh = opts->nhchainlength;
79 mstepsi = mstepsj = ns;
81 /* if scalefac is NULL, we are doing the NHC of the barostat */
84 if (scalefac == NULL) {
88 for (i=0; i<nvar; i++)
91 /* make it easier to iterate by selecting
92 out the sub-array that corresponds to this T group */
97 iQinv = &(MassQ->QPinv[i*nh]);
98 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
99 reft = max(0.0,opts->ref_t[0]);
100 Ekin = sqr(*veta)/MassQ->Winv;
102 iQinv = &(MassQ->Qinv[i*nh]);
103 tcstat = &ekind->tcstat[i];
105 reft = max(0.0,opts->ref_t[i]);
108 Ekin = 2*trace(tcstat->ekinf)*tcstat->ekinscalef_nhc;
110 Ekin = 2*trace(tcstat->ekinh)*tcstat->ekinscaleh_nhc;
115 for(mi=0;mi<mstepsi;mi++)
117 for(mj=0;mj<mstepsj;mj++)
119 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
120 dt = sy_const[ns][mj] * dtfull / mstepsi;
122 /* compute the thermal forces */
123 GQ[0] = iQinv[0]*(Ekin - nd*kT);
127 if (iQinv[j+1] > 0) {
128 /* we actually don't need to update here if we save the
129 state of the GQ, but it's easier to just recompute*/
130 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
136 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
139 Efac = exp(-0.125*dt*ivxi[j]);
140 ivxi[j-1] = Efac*(ivxi[j-1]*Efac + 0.25*dt*GQ[j-1]);
143 Efac = exp(-0.5*dt*ivxi[0]);
151 /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
152 able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
153 think about this a bit more . . . */
155 GQ[0] = iQinv[0]*(Ekin - nd*kT);
157 /* update thermostat positions */
160 ixi[j] += 0.5*dt*ivxi[j];
165 Efac = exp(-0.125*dt*ivxi[j+1]);
166 ivxi[j] = Efac*(ivxi[j]*Efac + 0.25*dt*GQ[j]);
167 if (iQinv[j+1] > 0) {
168 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
173 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
180 static void boxv_trotter(t_inputrec *ir, real *veta, real dt, tensor box,
181 gmx_ekindata_t *ekind, tensor vir, real pcorr, t_extmass *MassQ)
188 tensor Winvm,ekinmod,localpres;
190 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
191 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
194 if (ir->epct==epctSEMIISOTROPIC)
203 /* eta is in pure units. veta is in units of ps^-1. GW is in
204 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
205 taken to use only RATIOS of eta in updating the volume. */
207 /* we take the partial pressure tensors, modify the
208 kinetic energy tensor, and recovert to pressure */
210 if (ir->opts.nrdf[0]==0)
212 gmx_fatal(FARGS,"Barostat is coupled to a T-group with no degrees of freedom\n");
214 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
215 alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
216 alpha *= ekind->tcstat[0].ekinscalef_nhc;
217 msmul(ekind->ekin,alpha,ekinmod);
218 /* for now, we use Elr = 0, because if you want to get it right, you
219 really should be using PME. Maybe print a warning? */
221 pscal = calc_pres(ir->ePBC,nwall,box,ekinmod,vir,localpres)+pcorr;
224 GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
230 * This file implements temperature and pressure coupling algorithms:
231 * For now only the Weak coupling and the modified weak coupling.
233 * Furthermore computation of pressure and temperature is done here
237 real calc_pres(int ePBC,int nwall,matrix box,tensor ekin,tensor vir,
243 if (ePBC==epbcNONE || (ePBC==epbcXY && nwall!=2))
246 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
247 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
251 fac=PRESFAC*2.0/det(box);
252 for(n=0; (n<DIM); n++)
253 for(m=0; (m<DIM); m++)
254 pres[n][m] = (ekin[n][m] - vir[n][m])*fac;
257 pr_rvecs(debug,0,"PC: pres",pres,DIM);
258 pr_rvecs(debug,0,"PC: ekin",ekin,DIM);
259 pr_rvecs(debug,0,"PC: vir ",vir, DIM);
260 pr_rvecs(debug,0,"PC: box ",box, DIM);
263 return trace(pres)/DIM;
266 real calc_temp(real ekin,real nrdf)
269 return (2.0*ekin)/(nrdf*BOLTZ);
274 void parrinellorahman_pcoupl(FILE *fplog,gmx_large_int_t step,
275 t_inputrec *ir,real dt,tensor pres,
276 tensor box,tensor box_rel,tensor boxv,
277 tensor M,matrix mu,gmx_bool bFirstStep)
279 /* This doesn't do any coordinate updating. It just
280 * integrates the box vector equations from the calculated
281 * acceleration due to pressure difference. We also compute
282 * the tensor M which is used in update to couple the particle
283 * coordinates to the box vectors.
285 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
288 * M_nk = (h') * (h' * h + h' h) * h
290 * with the dots denoting time derivatives and h is the transformation from
291 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
292 * This also goes for the pressure and M tensors - they are transposed relative
293 * to ours. Our equation thus becomes:
296 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
298 * where b is the gromacs box matrix.
299 * Our box accelerations are given by
301 * b = vol/W inv(box') * (P-ref_P) (=h')
306 real vol=box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
307 real atot,arel,change,maxchange,xy_pressure;
308 tensor invbox,pdiff,t1,t2;
312 m_inv_ur0(box,invbox);
315 /* Note that PRESFAC does not occur here.
316 * The pressure and compressibility always occur as a product,
317 * therefore the pressure unit drops out.
319 maxl=max(box[XX][XX],box[YY][YY]);
320 maxl=max(maxl,box[ZZ][ZZ]);
324 (4*M_PI*M_PI*ir->compress[d][n])/(3*ir->tau_p*ir->tau_p*maxl);
326 m_sub(pres,ir->ref_p,pdiff);
328 if(ir->epct==epctSURFACETENSION) {
329 /* Unlike Berendsen coupling it might not be trivial to include a z
330 * pressure correction here? On the other hand we don't scale the
331 * box momentarily, but change accelerations, so it might not be crucial.
333 xy_pressure=0.5*(pres[XX][XX]+pres[YY][YY]);
335 pdiff[d][d]=(xy_pressure-(pres[ZZ][ZZ]-ir->ref_p[d][d]/box[d][d]));
338 tmmul(invbox,pdiff,t1);
339 /* Move the off-diagonal elements of the 'force' to one side to ensure
340 * that we obey the box constraints.
344 t1[d][n] += t1[n][d];
350 case epctANISOTROPIC:
353 t1[d][n] *= winv[d][n]*vol;
356 /* calculate total volume acceleration */
357 atot=box[XX][XX]*box[YY][YY]*t1[ZZ][ZZ]+
358 box[XX][XX]*t1[YY][YY]*box[ZZ][ZZ]+
359 t1[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
361 /* set all RELATIVE box accelerations equal, and maintain total V
365 t1[d][n] = winv[0][0]*vol*arel*box[d][n];
367 case epctSEMIISOTROPIC:
368 case epctSURFACETENSION:
369 /* Note the correction to pdiff above for surftens. coupling */
371 /* calculate total XY volume acceleration */
372 atot=box[XX][XX]*t1[YY][YY]+t1[XX][XX]*box[YY][YY];
373 arel=atot/(2*box[XX][XX]*box[YY][YY]);
374 /* set RELATIVE XY box accelerations equal, and maintain total V
375 * change speed. Dont change the third box vector accelerations */
378 t1[d][n] = winv[d][n]*vol*arel*box[d][n];
380 t1[ZZ][n] *= winv[d][n]*vol;
383 gmx_fatal(FARGS,"Parrinello-Rahman pressure coupling type %s "
384 "not supported yet\n",EPCOUPLTYPETYPE(ir->epct));
391 boxv[d][n] += dt*t1[d][n];
393 /* We do NOT update the box vectors themselves here, since
394 * we need them for shifting later. It is instead done last
395 * in the update() routine.
398 /* Calculate the change relative to diagonal elements-
399 since it's perfectly ok for the off-diagonal ones to
400 be zero it doesn't make sense to check the change relative
404 change=fabs(dt*boxv[d][n]/box[d][d]);
406 if (change>maxchange)
410 if (maxchange > 0.01 && fplog) {
413 "\nStep %s Warning: Pressure scaling more than 1%%. "
414 "This may mean your system\n is not yet equilibrated. "
415 "Use of Parrinello-Rahman pressure coupling during\n"
416 "equilibration can lead to simulation instability, "
417 "and is discouraged.\n",
418 gmx_step_str(step,buf));
422 preserve_box_shape(ir,box_rel,boxv);
424 mtmul(boxv,box,t1); /* t1=boxv * b' */
428 /* Determine the scaling matrix mu for the coordinates */
431 t1[d][n] = box[d][n] + dt*boxv[d][n];
432 preserve_box_shape(ir,box_rel,t1);
433 /* t1 is the box at t+dt, determine mu as the relative change */
434 mmul_ur0(invbox,t1,mu);
437 void berendsen_pcoupl(FILE *fplog,gmx_large_int_t step,
438 t_inputrec *ir,real dt, tensor pres,matrix box,
442 real scalar_pressure, xy_pressure, p_corr_z;
443 char *ptr,buf[STRLEN];
446 * Calculate the scaling matrix mu
450 for(d=0; d<DIM; d++) {
451 scalar_pressure += pres[d][d]/DIM;
453 xy_pressure += pres[d][d]/(DIM-1);
455 /* Pressure is now in bar, everywhere. */
456 #define factor(d,m) (ir->compress[d][m]*dt/ir->tau_p)
458 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
459 * necessary for triclinic scaling
466 mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
469 case epctSEMIISOTROPIC:
471 mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
473 1.0 - factor(ZZ,ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
475 case epctANISOTROPIC:
478 mu[d][n] = (d==n ? 1.0 : 0.0)
479 -factor(d,n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
481 case epctSURFACETENSION:
482 /* ir->ref_p[0/1] is the reference surface-tension times *
483 * the number of surfaces */
484 if (ir->compress[ZZ][ZZ])
485 p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
487 /* when the compressibity is zero, set the pressure correction *
488 * in the z-direction to zero to get the correct surface tension */
490 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
491 for(d=0; d<DIM-1; d++)
492 mu[d][d] = 1.0 + factor(d,d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
493 - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
496 gmx_fatal(FARGS,"Berendsen pressure coupling type %s not supported yet\n",
497 EPCOUPLTYPETYPE(ir->epct));
500 /* To fullfill the orientation restrictions on triclinic boxes
501 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
502 * the other elements of mu to first order.
504 mu[YY][XX] += mu[XX][YY];
505 mu[ZZ][XX] += mu[XX][ZZ];
506 mu[ZZ][YY] += mu[YY][ZZ];
512 pr_rvecs(debug,0,"PC: pres ",pres,3);
513 pr_rvecs(debug,0,"PC: mu ",mu,3);
516 if (mu[XX][XX]<0.99 || mu[XX][XX]>1.01 ||
517 mu[YY][YY]<0.99 || mu[YY][YY]>1.01 ||
518 mu[ZZ][ZZ]<0.99 || mu[ZZ][ZZ]>1.01) {
520 sprintf(buf,"\nStep %s Warning: pressure scaling more than 1%%, "
522 gmx_step_str(step,buf2),mu[XX][XX],mu[YY][YY],mu[ZZ][ZZ]);
524 fprintf(fplog,"%s",buf);
525 fprintf(stderr,"%s",buf);
529 void berendsen_pscale(t_inputrec *ir,matrix mu,
530 matrix box,matrix box_rel,
531 int start,int nr_atoms,
532 rvec x[],unsigned short cFREEZE[],
535 ivec *nFreeze=ir->opts.nFreeze;
538 /* Scale the positions */
539 for (n=start; n<start+nr_atoms; n++) {
544 x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
546 x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
548 x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
550 /* compute final boxlengths */
551 for (d=0; d<DIM; d++) {
552 box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
553 box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
554 box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
557 preserve_box_shape(ir,box_rel,box);
559 /* (un)shifting should NOT be done after this,
560 * since the box vectors might have changed
562 inc_nrnb(nrnb,eNR_PCOUPL,nr_atoms);
565 void berendsen_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt)
573 for(i=0; (i<opts->ngtc); i++)
577 T = ekind->tcstat[i].T;
581 T = ekind->tcstat[i].Th;
584 if ((opts->tau_t[i] > 0) && (T > 0.0)) {
585 reft = max(0.0,opts->ref_t[i]);
586 lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
587 ekind->tcstat[i].lambda = max(min(lll,1.25),0.8);
590 ekind->tcstat[i].lambda = 1.0;
595 fprintf(debug,"TC: group %d: T: %g, Lambda: %g\n",
596 i,T,ekind->tcstat[i].lambda);
601 static int poisson_variate(real lambda,gmx_rng_t rng) {
612 p *= gmx_rng_uniform_real(rng);
618 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)
621 int i,j,k,d,len,n,ngtc,gc=0;
622 int nshake, nsettle, nrandom, nrand_group;
623 real boltz,scal,reft,prand;
626 /* convenience variables */
630 /* idef is only passed in if it's chance-per-particle andersen, so
631 it essentially serves as a boolean to determine which type of
632 andersen is being used */
635 /* randomly atoms to randomize. However, all constraint
636 groups have to have either all of the atoms or none of the
640 1. Select whether or not to randomize each atom to get the correct probability.
641 2. Cycle through the constraint groups.
642 2a. for each constraint group, determine the fraction f of that constraint group that are
643 chosen to be randomized.
644 2b. all atoms in the constraint group are randomized with probability f.
648 if ((rate < 0.05) && (md->homenr > 50))
650 /* if the rate is relatively high, use a standard method, if low rate,
652 /* poisson distributions approxmation, more efficient for
653 * low rates, fewer random numbers required */
654 nrandom = poisson_variate(md->homenr*rate,rng); /* how many do we randomize? Use Poisson. */
655 /* now we know how many, choose them randomly. No worries about repeats, at this rate, it's negligible.
656 worst thing that happens, it lowers the true rate an negligible amount */
657 for (i=0;i<nrandom;i++)
659 randatom[(int)(gmx_rng_uniform_real(rng)*md->homenr)] = TRUE;
664 for (i=0;i<md->homenr;i++)
666 if (gmx_rng_uniform_real(rng)<rate)
674 /* instead of looping over the constraint groups, if we had a
675 list of which atoms were in which constraint groups, we
676 could then loop over only the groups that are randomized
677 now. But that is not available now. Create later after
678 determining whether there actually is any slowing. */
680 /* first, loop through the settles to make sure all groups either entirely randomized, or not randomized. */
682 nsettle = idef->il[F_SETTLE].nr/2;
683 for (i=0;i<nsettle;i++)
685 iatoms = idef->il[F_SETTLE].iatoms;
687 for (k=0;k<3;k++) /* settles are always 3 atoms, hardcoded */
689 if (randatom[iatoms[2*i+1]+k])
691 nrand_group++; /* count the number of atoms to be shaken in the settles group */
692 randatom[iatoms[2*i+1]+k] = FALSE;
698 prand = (nrand_group)/3.0; /* use this fraction to compute the probability the
699 whole group is randomized */
700 if (gmx_rng_uniform_real(rng)<prand)
704 randatom[iatoms[2*i+1]+k] = TRUE; /* mark them all to be randomized */
711 /* now loop through the shake groups */
713 for (i=0;i<nshake;i++)
715 iatoms = &(idef->il[F_CONSTR].iatoms[sblock[i]]);
716 len = sblock[i+1]-sblock[i];
721 { /* only 2/3 of the sblock items are atoms, the others are labels */
722 if (randatom[iatoms[k]])
725 randatom[iatoms[k]] = FALSE; /* need to mark it false here in case the atom is in more than
726 one group in the shake block */
733 prand = (nrand_group)/(1.0*(2*len/3));
734 if (gmx_rng_uniform_real(rng)<prand)
739 { /* only 2/3 of the sblock items are atoms, the others are labels */
740 randatom[iatoms[k]] = TRUE; /* randomize all of them */
750 for (i=0;i<md->homenr;i++) /* now loop over the list of atoms */
754 randatom_list[n] = i;
758 nrandom = n; /* there are some values of nrandom for which
759 this algorithm won't work; for example all
760 water molecules and nrandom =/= 3. Better to
761 recount and use this number (which we
762 calculate anyway: it will not affect
763 the average number of atoms accepted.
769 /* if it's andersen-massive, then randomize all the atoms */
770 nrandom = md->homenr;
771 for (i=0;i<nrandom;i++)
773 randatom_list[i] = i;
777 /* randomize the velocities of the selected particles */
779 for (i=0;i<nrandom;i++) /* now loop over the list of atoms */
781 n = randatom_list[i];
784 gc = md->cTC[n]; /* assign the atom to a temperature group if there are more than one */
788 scal = sqrt(boltzfac[gc]*md->invmass[n]);
791 state->v[n][d] = scal*gmx_rng_gaussian_table(rng);
794 randatom[n] = FALSE; /* unmark this atom for randomization */
799 void nosehoover_tcoupl(t_grpopts *opts,gmx_ekindata_t *ekind,real dt,
800 double xi[],double vxi[], t_extmass *MassQ)
805 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
807 for(i=0; (i<opts->ngtc); i++)
809 reft = max(0.0,opts->ref_t[i]);
811 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
812 xi[i] += dt*(oldvxi + vxi[i])*0.5;
816 t_state *init_bufstate(const t_state *template_state)
819 int nc = template_state->nhchainlength;
821 snew(state->nosehoover_xi,nc*template_state->ngtc);
822 snew(state->nosehoover_vxi,nc*template_state->ngtc);
823 snew(state->therm_integral,template_state->ngtc);
824 snew(state->nhpres_xi,nc*template_state->nnhpres);
825 snew(state->nhpres_vxi,nc*template_state->nnhpres);
830 void destroy_bufstate(t_state *state)
834 sfree(state->nosehoover_xi);
835 sfree(state->nosehoover_vxi);
836 sfree(state->therm_integral);
837 sfree(state->nhpres_xi);
838 sfree(state->nhpres_vxi);
842 void trotter_update(t_inputrec *ir,gmx_large_int_t step, gmx_ekindata_t *ekind,
843 gmx_enerdata_t *enerd, t_state *state,
844 tensor vir, t_mdatoms *md,
845 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
848 int n,i,j,d,ntgrp,ngtc,gc=0;
849 t_grp_tcstat *tcstat;
851 gmx_large_int_t step_eff;
852 real ecorr,pcorr,dvdlcorr;
853 real bmass,qmass,reft,kT,dt,nd;
854 tensor dumpres,dumvir;
855 double *scalefac,dtc;
857 rvec sumv={0,0,0},consk;
860 if (trotter_seqno <= ettTSEQ2)
862 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
863 is actually the last half step from the previous step. Thus the first half step
864 actually corresponds to the n-1 step*/
870 bCouple = (ir->nsttcouple == 1 ||
871 do_per_step(step_eff+ir->nsttcouple,ir->nsttcouple));
873 trotter_seq = trotter_seqlist[trotter_seqno];
875 if ((trotter_seq[0] == etrtSKIPALL) || (!bCouple))
879 dtc = ir->nsttcouple*ir->delta_t; /* This is OK for NPT, because nsttcouple == nstpcouple is enforcesd */
880 opts = &(ir->opts); /* just for ease of referencing */
883 snew(scalefac,opts->ngtc);
888 /* execute the series of trotter updates specified in the trotterpart array */
890 for (i=0;i<NTROTTERPARTS;i++){
891 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
892 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
901 switch (trotter_seq[i])
905 boxv_trotter(ir,&(state->veta),dt,state->box,ekind,vir,
906 enerd->term[F_PDISPCORR],MassQ);
910 NHC_trotter(opts,state->nnhpres,ekind,dt,state->nhpres_xi,
911 state->nhpres_vxi,NULL,&(state->veta),MassQ,FALSE);
915 NHC_trotter(opts,opts->ngtc,ekind,dt,state->nosehoover_xi,
916 state->nosehoover_vxi,scalefac,NULL,MassQ,(ir->eI==eiVV));
917 /* need to rescale the kinetic energies and velocities here. Could
918 scale the velocities later, but we need them scaled in order to
919 produce the correct outputs, so we'll scale them here. */
921 for (i=0; i<ngtc;i++)
923 tcstat = &ekind->tcstat[i];
924 tcstat->vscale_nhc = scalefac[i];
925 tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
926 tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
928 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
929 /* but do we actually need the total? */
931 /* modify the velocities as well */
932 for (n=md->start;n<md->start+md->homenr;n++)
934 if (md->cTC) /* does this conditional need to be here? is this always true?*/
940 state->v[n][d] *= scalefac[gc];
947 sumv[d] += (state->v[n][d])/md->invmass[n];
956 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
964 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
966 fprintf(debug,"Conserved kappa: %15.8f %15.8f %15.8f\n",consk[0],consk[1],consk[2]);
974 extern void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit)
976 int n,i,j,d,ntgrp,ngtc,nnhpres,nh,gc=0;
977 t_grp_tcstat *tcstat;
979 real ecorr,pcorr,dvdlcorr;
980 real bmass,qmass,reft,kT,dt,ndj,nd;
981 tensor dumpres,dumvir;
983 opts = &(ir->opts); /* just for ease of referencing */
984 ngtc = ir->opts.ngtc;
985 nnhpres = state->nnhpres;
986 nh = state->nhchainlength;
988 if (ir->eI == eiMD) {
991 snew(MassQ->Qinv,ngtc);
993 for(i=0; (i<ngtc); i++)
995 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
997 MassQ->Qinv[i]=1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
1005 else if (EI_VV(ir->eI))
1007 /* Set pressure variables */
1011 if (state->vol0 == 0)
1013 state->vol0 = det(state->box);
1014 /* because we start by defining a fixed
1015 compressibility, we need the volume at this
1016 compressibility to solve the problem. */
1020 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1021 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1022 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
1023 /* An alternate mass definition, from Tuckerman et al. */
1024 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1029 MassQ->Winvm[d][n]= PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
1030 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1031 before using MTTK for anisotropic states.*/
1034 /* Allocate space for thermostat variables */
1037 snew(MassQ->Qinv,ngtc*nh);
1040 /* now, set temperature variables */
1041 for (i=0; i<ngtc; i++)
1043 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1045 reft = max(0.0,opts->ref_t[i]);
1058 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
1066 MassQ->Qinv[i*nh+j] = 0.0;
1073 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
1075 int n,i,j,d,ntgrp,ngtc,nnhpres,nh,gc=0;
1076 t_grp_tcstat *tcstat;
1078 real ecorr,pcorr,dvdlcorr;
1079 real bmass,qmass,reft,kT,dt,ndj,nd;
1080 tensor dumpres,dumvir;
1083 opts = &(ir->opts); /* just for ease of referencing */
1085 nnhpres = state->nnhpres;
1086 nh = state->nhchainlength;
1088 init_npt_masses(ir, state, MassQ, TRUE);
1090 /* first, initialize clear all the trotter calls */
1091 snew(trotter_seq,ettTSEQMAX);
1092 for (i=0;i<ettTSEQMAX;i++)
1094 snew(trotter_seq[i],NTROTTERPARTS);
1095 for (j=0;j<NTROTTERPARTS;j++) {
1096 trotter_seq[i][j] = etrtNONE;
1098 trotter_seq[i][0] = etrtSKIPALL;
1103 /* no trotter calls, so we never use the values in the array.
1104 * We access them (so we need to define them, but ignore
1110 /* compute the kinetic energy by using the half step velocities or
1111 * the kinetic energies, depending on the order of the trotter calls */
1115 if (IR_NPT_TROTTER(ir))
1117 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1118 We start with the initial one. */
1119 /* first, a round that estimates veta. */
1120 trotter_seq[0][0] = etrtBAROV;
1122 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1124 /* The first half trotter update */
1125 trotter_seq[2][0] = etrtBAROV;
1126 trotter_seq[2][1] = etrtNHC;
1127 trotter_seq[2][2] = etrtBARONHC;
1129 /* The second half trotter update */
1130 trotter_seq[3][0] = etrtBARONHC;
1131 trotter_seq[3][1] = etrtNHC;
1132 trotter_seq[3][2] = etrtBAROV;
1134 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1137 else if (IR_NVT_TROTTER(ir))
1139 /* This is the easy version - there are only two calls, both the same.
1140 Otherwise, even easier -- no calls */
1141 trotter_seq[2][0] = etrtNHC;
1142 trotter_seq[3][0] = etrtNHC;
1144 else if (IR_NPH_TROTTER(ir))
1146 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1147 We start with the initial one. */
1148 /* first, a round that estimates veta. */
1149 trotter_seq[0][0] = etrtBAROV;
1151 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1153 /* The first half trotter update */
1154 trotter_seq[2][0] = etrtBAROV;
1155 trotter_seq[2][1] = etrtBARONHC;
1157 /* The second half trotter update */
1158 trotter_seq[3][0] = etrtBARONHC;
1159 trotter_seq[3][1] = etrtBAROV;
1161 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1164 else if (ir->eI==eiVVAK)
1166 if (IR_NPT_TROTTER(ir))
1168 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1169 We start with the initial one. */
1170 /* first, a round that estimates veta. */
1171 trotter_seq[0][0] = etrtBAROV;
1173 /* The first half trotter update, part 1 -- double update, because it commutes */
1174 trotter_seq[1][0] = etrtNHC;
1176 /* The first half trotter update, part 2 */
1177 trotter_seq[2][0] = etrtBAROV;
1178 trotter_seq[2][1] = etrtBARONHC;
1180 /* The second half trotter update, part 1 */
1181 trotter_seq[3][0] = etrtBARONHC;
1182 trotter_seq[3][1] = etrtBAROV;
1184 /* The second half trotter update */
1185 trotter_seq[4][0] = etrtNHC;
1187 else if (IR_NVT_TROTTER(ir))
1189 /* This is the easy version - there is only one call, both the same.
1190 Otherwise, even easier -- no calls */
1191 trotter_seq[1][0] = etrtNHC;
1192 trotter_seq[4][0] = etrtNHC;
1194 else if (IR_NPH_TROTTER(ir))
1196 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1197 We start with the initial one. */
1198 /* first, a round that estimates veta. */
1199 trotter_seq[0][0] = etrtBAROV;
1201 /* The first half trotter update, part 1 -- leave zero */
1202 trotter_seq[1][0] = etrtNHC;
1204 /* The first half trotter update, part 2 */
1205 trotter_seq[2][0] = etrtBAROV;
1206 trotter_seq[2][1] = etrtBARONHC;
1208 /* The second half trotter update, part 1 */
1209 trotter_seq[3][0] = etrtBARONHC;
1210 trotter_seq[3][1] = etrtBAROV;
1212 /* The second half trotter update -- blank for now */
1220 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
1223 snew(MassQ->QPinv,nnhpres*opts->nhchainlength);
1225 /* barostat temperature */
1226 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1228 reft = max(0.0,opts->ref_t[0]);
1230 for (i=0;i<nnhpres;i++) {
1240 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1246 for (i=0;i<nnhpres;i++) {
1249 MassQ->QPinv[i*nh+j]=0.0;
1256 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1258 int i,j,nd,ndj,bmass,qmass,ngtcall;
1259 real ener_npt,reft,eta,kT,tau;
1263 int nh = state->nhchainlength;
1267 /* now we compute the contribution of the pressure to the conserved quantity*/
1269 if (ir->epc==epcMTTK)
1271 /* find the volume, and the kinetic energy of the volume */
1276 /* contribution from the pressure momenenta */
1277 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1279 /* contribution from the PV term */
1280 vol = det(state->box);
1281 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1284 case epctANISOTROPIC:
1288 case epctSURFACETENSION:
1291 case epctSEMIISOTROPIC:
1299 if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1301 /* add the energy from the barostat thermostat chain */
1302 for (i=0;i<state->nnhpres;i++) {
1304 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1305 ivxi = &state->nhpres_vxi[i*nh];
1306 ixi = &state->nhpres_xi[i*nh];
1307 iQinv = &(MassQ->QPinv[i*nh]);
1308 reft = max(ir->opts.ref_t[0],0); /* using 'System' temperature */
1315 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1316 /* contribution from the thermal variable of the NH chain */
1317 ener_npt += ixi[j]*kT;
1321 fprintf(debug,"P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f",i,j,ivxi[j],ixi[j]);
1329 for(i=0; i<ir->opts.ngtc; i++)
1331 ixi = &state->nosehoover_xi[i*nh];
1332 ivxi = &state->nosehoover_vxi[i*nh];
1333 iQinv = &(MassQ->Qinv[i*nh]);
1335 nd = ir->opts.nrdf[i];
1336 reft = max(ir->opts.ref_t[i],0);
1341 if (IR_NVT_TROTTER(ir))
1343 /* contribution from the thermal momenta of the NH chain */
1347 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1348 /* contribution from the thermal variable of the NH chain */
1356 ener_npt += ndj*ixi[j]*kT;
1360 else /* Other non Trotter temperature NH control -- no chains yet. */
1362 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1363 ener_npt += nd*ixi[0]*kT;
1371 static real vrescale_gamdev(int ia, gmx_rng_t rng)
1372 /* Gamma distribution, adapted from numerical recipes */
1375 real am,e,s,v1,v2,x,y;
1382 for(j=1; j<=ia; j++)
1384 x *= gmx_rng_uniform_real(rng);
1398 v1 = gmx_rng_uniform_real(rng);
1399 v2 = 2.0*gmx_rng_uniform_real(rng)-1.0;
1401 while (v1*v1 + v2*v2 > 1.0 ||
1402 v1*v1*GMX_REAL_MAX < 3.0*ia);
1403 /* The last check above ensures that both x (3.0 > 2.0 in s)
1404 * and the pre-factor for e do not go out of range.
1408 s = sqrt(2.0*am + 1.0);
1412 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1414 while (gmx_rng_uniform_real(rng) > e);
1420 static real vrescale_sumnoises(int nn,gmx_rng_t rng)
1423 * Returns the sum of n independent gaussian noises squared
1424 * (i.e. equivalent to summing the square of the return values
1425 * of nn calls to gmx_rng_gaussian_real).xs
1431 } else if (nn == 1) {
1432 rr = gmx_rng_gaussian_real(rng);
1434 } else if (nn % 2 == 0) {
1435 return 2.0*vrescale_gamdev(nn/2,rng);
1437 rr = gmx_rng_gaussian_real(rng);
1438 return 2.0*vrescale_gamdev((nn-1)/2,rng) + rr*rr;
1442 static real vrescale_resamplekin(real kk,real sigma, int ndeg, real taut,
1446 * Generates a new value for the kinetic energy,
1447 * according to Bussi et al JCP (2007), Eq. (A7)
1448 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1449 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1450 * ndeg: number of degrees of freedom of the atoms to be thermalized
1451 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1456 factor = exp(-1.0/taut);
1460 rr = gmx_rng_gaussian_real(rng);
1463 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1,rng) + rr*rr)/ndeg - kk) +
1464 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1467 void vrescale_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt,
1468 double therm_integral[],gmx_rng_t rng)
1472 real Ek,Ek_ref1,Ek_ref,Ek_new;
1476 for(i=0; (i<opts->ngtc); i++)
1480 Ek = trace(ekind->tcstat[i].ekinf);
1484 Ek = trace(ekind->tcstat[i].ekinh);
1487 if (opts->tau_t[i] > 0 && opts->nrdf[i] > 0 && Ek > 0)
1489 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1490 Ek_ref = Ek_ref1*opts->nrdf[i];
1492 Ek_new = vrescale_resamplekin(Ek,Ek_ref,opts->nrdf[i],
1493 opts->tau_t[i]/dt,rng);
1495 /* Analytically Ek_new>=0, but we check for rounding errors */
1498 ekind->tcstat[i].lambda = 0.0;
1502 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1505 therm_integral[i] -= Ek_new - Ek;
1509 fprintf(debug,"TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1510 i,Ek_ref,Ek,Ek_new,ekind->tcstat[i].lambda);
1515 ekind->tcstat[i].lambda = 1.0;
1520 real vrescale_energy(t_grpopts *opts,double therm_integral[])
1526 for(i=0; i<opts->ngtc; i++) {
1527 ener += therm_integral[i];
1533 void rescale_velocities(gmx_ekindata_t *ekind,t_mdatoms *mdatoms,
1534 int start,int end,rvec v[])
1537 t_grp_tcstat *tcstat;
1538 unsigned short *cACC,*cTC;
1543 tcstat = ekind->tcstat;
1548 gstat = ekind->grpstat;
1549 cACC = mdatoms->cACC;
1553 for(n=start; n<end; n++)
1563 /* Only scale the velocity component relative to the COM velocity */
1564 rvec_sub(v[n],gstat[ga].u,vrel);
1565 lg = tcstat[gt].lambda;
1566 for(d=0; d<DIM; d++)
1568 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1575 for(n=start; n<end; n++)
1581 lg = tcstat[gt].lambda;
1582 for(d=0; d<DIM; d++)
1591 /* set target temperatures if we are annealing */
1592 void update_annealing_target_temp(t_grpopts *opts,real t)
1595 real pert,thist=0,x;
1597 for(i=0;i<opts->ngtc;i++) {
1598 npoints = opts->anneal_npoints[i];
1599 switch (opts->annealing[i]) {
1603 /* calculate time modulo the period */
1604 pert = opts->anneal_time[i][npoints-1];
1606 thist = t - n*pert; /* modulo time */
1607 /* Make sure rounding didn't get us outside the interval */
1608 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1615 gmx_fatal(FARGS,"Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",i,opts->ngtc,npoints);
1617 /* We are doing annealing for this group if we got here,
1618 * and we have the (relative) time as thist.
1619 * calculate target temp */
1621 while ((j < npoints-1) && (thist>(opts->anneal_time[i][j+1])))
1623 if (j < npoints-1) {
1624 /* Found our position between points j and j+1.
1625 * Interpolate: x is the amount from j+1, (1-x) from point j
1626 * First treat possible jumps in temperature as a special case.
1628 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1629 opts->ref_t[i]=opts->anneal_temp[i][j+1];
1631 x = ((thist-opts->anneal_time[i][j])/
1632 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1633 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1637 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];