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 /* these integration routines are only referenced inside this file */
57 static void NHC_trotter(t_grpopts *opts,int nvar, gmx_ekindata_t *ekind,real dtfull,
58 double xi[],double vxi[], double scalefac[], real *veta, t_extmass *MassQ, gmx_bool bEkinAveVel)
61 /* general routine for both barostat and thermostat nose hoover chains */
64 double Ekin,Efac,reft,kT,nd;
72 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
73 int nh = opts->nhchainlength;
76 mstepsi = mstepsj = ns;
78 /* if scalefac is NULL, we are doing the NHC of the barostat */
81 if (scalefac == NULL) {
85 for (i=0; i<nvar; i++)
88 /* make it easier to iterate by selecting
89 out the sub-array that corresponds to this T group */
94 iQinv = &(MassQ->QPinv[i*nh]);
95 nd = 1.0; /* THIS WILL CHANGE IF NOT ISOTROPIC */
96 reft = max(0.0,opts->ref_t[0]);
97 Ekin = sqr(*veta)/MassQ->Winv;
99 iQinv = &(MassQ->Qinv[i*nh]);
100 tcstat = &ekind->tcstat[i];
102 reft = max(0.0,opts->ref_t[i]);
105 Ekin = 2*trace(tcstat->ekinf)*tcstat->ekinscalef_nhc;
107 Ekin = 2*trace(tcstat->ekinh)*tcstat->ekinscaleh_nhc;
112 for(mi=0;mi<mstepsi;mi++)
114 for(mj=0;mj<mstepsj;mj++)
116 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
117 dt = sy_const[ns][mj] * dtfull / mstepsi;
119 /* compute the thermal forces */
120 GQ[0] = iQinv[0]*(Ekin - nd*kT);
124 if (iQinv[j+1] > 0) {
125 /* we actually don't need to update here if we save the
126 state of the GQ, but it's easier to just recompute*/
127 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
133 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
136 Efac = exp(-0.125*dt*ivxi[j]);
137 ivxi[j-1] = Efac*(ivxi[j-1]*Efac + 0.25*dt*GQ[j-1]);
140 Efac = exp(-0.5*dt*ivxi[0]);
148 /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
149 able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
150 think about this a bit more . . . */
152 GQ[0] = iQinv[0]*(Ekin - nd*kT);
154 /* update thermostat positions */
157 ixi[j] += 0.5*dt*ivxi[j];
162 Efac = exp(-0.125*dt*ivxi[j+1]);
163 ivxi[j] = Efac*(ivxi[j]*Efac + 0.25*dt*GQ[j]);
164 if (iQinv[j+1] > 0) {
165 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
170 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
177 static void boxv_trotter(t_inputrec *ir, real *veta, real dt, tensor box,
178 gmx_ekindata_t *ekind, tensor vir, real pcorr, real ecorr, t_extmass *MassQ)
185 tensor Winvm,ekinmod,localpres;
187 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
188 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
191 if (ir->epct==epctSEMIISOTROPIC)
200 /* eta is in pure units. veta is in units of ps^-1. GW is in
201 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
202 taken to use only RATIOS of eta in updating the volume. */
204 /* we take the partial pressure tensors, modify the
205 kinetic energy tensor, and recovert to pressure */
207 if (ir->opts.nrdf[0]==0)
209 gmx_fatal(FARGS,"Barostat is coupled to a T-group with no degrees of freedom\n");
211 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
212 alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
213 alpha *= ekind->tcstat[0].ekinscalef_nhc;
214 msmul(ekind->ekin,alpha,ekinmod);
215 /* for now, we use Elr = 0, because if you want to get it right, you
216 really should be using PME. Maybe print a warning? */
218 pscal = calc_pres(ir->ePBC,nwall,box,ekinmod,vir,localpres);
221 GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
227 * This file implements temperature and pressure coupling algorithms:
228 * For now only the Weak coupling and the modified weak coupling.
230 * Furthermore computation of pressure and temperature is done here
234 real calc_pres(int ePBC,int nwall,matrix box,tensor ekin,tensor vir,
240 if (ePBC==epbcNONE || (ePBC==epbcXY && nwall!=2))
243 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
244 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
248 fac=PRESFAC*2.0/det(box);
249 for(n=0; (n<DIM); n++)
250 for(m=0; (m<DIM); m++)
251 pres[n][m] = (ekin[n][m] - vir[n][m])*fac;
254 pr_rvecs(debug,0,"PC: pres",pres,DIM);
255 pr_rvecs(debug,0,"PC: ekin",ekin,DIM);
256 pr_rvecs(debug,0,"PC: vir ",vir, DIM);
257 pr_rvecs(debug,0,"PC: box ",box, DIM);
260 return trace(pres)/DIM;
263 real calc_temp(real ekin,real nrdf)
266 return (2.0*ekin)/(nrdf*BOLTZ);
271 void parrinellorahman_pcoupl(FILE *fplog,gmx_large_int_t step,
272 t_inputrec *ir,real dt,tensor pres,
273 tensor box,tensor box_rel,tensor boxv,
274 tensor M,matrix mu,gmx_bool bFirstStep)
276 /* This doesn't do any coordinate updating. It just
277 * integrates the box vector equations from the calculated
278 * acceleration due to pressure difference. We also compute
279 * the tensor M which is used in update to couple the particle
280 * coordinates to the box vectors.
282 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
285 * M_nk = (h') * (h' * h + h' h) * h
287 * with the dots denoting time derivatives and h is the transformation from
288 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
289 * This also goes for the pressure and M tensors - they are transposed relative
290 * to ours. Our equation thus becomes:
293 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
295 * where b is the gromacs box matrix.
296 * Our box accelerations are given by
298 * b = vol/W inv(box') * (P-ref_P) (=h')
303 real vol=box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
304 real atot,arel,change,maxchange,xy_pressure;
305 tensor invbox,pdiff,t1,t2;
309 m_inv_ur0(box,invbox);
312 /* Note that PRESFAC does not occur here.
313 * The pressure and compressibility always occur as a product,
314 * therefore the pressure unit drops out.
316 maxl=max(box[XX][XX],box[YY][YY]);
317 maxl=max(maxl,box[ZZ][ZZ]);
321 (4*M_PI*M_PI*ir->compress[d][n])/(3*ir->tau_p*ir->tau_p*maxl);
323 m_sub(pres,ir->ref_p,pdiff);
325 if(ir->epct==epctSURFACETENSION) {
326 /* Unlike Berendsen coupling it might not be trivial to include a z
327 * pressure correction here? On the other hand we don't scale the
328 * box momentarily, but change accelerations, so it might not be crucial.
330 xy_pressure=0.5*(pres[XX][XX]+pres[YY][YY]);
332 pdiff[d][d]=(xy_pressure-(pres[ZZ][ZZ]-ir->ref_p[d][d]/box[d][d]));
335 tmmul(invbox,pdiff,t1);
336 /* Move the off-diagonal elements of the 'force' to one side to ensure
337 * that we obey the box constraints.
341 t1[d][n] += t1[n][d];
347 case epctANISOTROPIC:
350 t1[d][n] *= winv[d][n]*vol;
353 /* calculate total volume acceleration */
354 atot=box[XX][XX]*box[YY][YY]*t1[ZZ][ZZ]+
355 box[XX][XX]*t1[YY][YY]*box[ZZ][ZZ]+
356 t1[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
358 /* set all RELATIVE box accelerations equal, and maintain total V
362 t1[d][n] = winv[0][0]*vol*arel*box[d][n];
364 case epctSEMIISOTROPIC:
365 case epctSURFACETENSION:
366 /* Note the correction to pdiff above for surftens. coupling */
368 /* calculate total XY volume acceleration */
369 atot=box[XX][XX]*t1[YY][YY]+t1[XX][XX]*box[YY][YY];
370 arel=atot/(2*box[XX][XX]*box[YY][YY]);
371 /* set RELATIVE XY box accelerations equal, and maintain total V
372 * change speed. Dont change the third box vector accelerations */
375 t1[d][n] = winv[d][n]*vol*arel*box[d][n];
377 t1[ZZ][n] *= winv[d][n]*vol;
380 gmx_fatal(FARGS,"Parrinello-Rahman pressure coupling type %s "
381 "not supported yet\n",EPCOUPLTYPETYPE(ir->epct));
388 boxv[d][n] += dt*t1[d][n];
390 /* We do NOT update the box vectors themselves here, since
391 * we need them for shifting later. It is instead done last
392 * in the update() routine.
395 /* Calculate the change relative to diagonal elements-
396 since it's perfectly ok for the off-diagonal ones to
397 be zero it doesn't make sense to check the change relative
401 change=fabs(dt*boxv[d][n]/box[d][d]);
403 if (change>maxchange)
407 if (maxchange > 0.01 && fplog) {
410 "\nStep %s Warning: Pressure scaling more than 1%%. "
411 "This may mean your system\n is not yet equilibrated. "
412 "Use of Parrinello-Rahman pressure coupling during\n"
413 "equilibration can lead to simulation instability, "
414 "and is discouraged.\n",
415 gmx_step_str(step,buf));
419 preserve_box_shape(ir,box_rel,boxv);
421 mtmul(boxv,box,t1); /* t1=boxv * b' */
425 /* Determine the scaling matrix mu for the coordinates */
428 t1[d][n] = box[d][n] + dt*boxv[d][n];
429 preserve_box_shape(ir,box_rel,t1);
430 /* t1 is the box at t+dt, determine mu as the relative change */
431 mmul_ur0(invbox,t1,mu);
434 void berendsen_pcoupl(FILE *fplog,gmx_large_int_t step,
435 t_inputrec *ir,real dt, tensor pres,matrix box,
439 real scalar_pressure, xy_pressure, p_corr_z;
440 char *ptr,buf[STRLEN];
443 * Calculate the scaling matrix mu
447 for(d=0; d<DIM; d++) {
448 scalar_pressure += pres[d][d]/DIM;
450 xy_pressure += pres[d][d]/(DIM-1);
452 /* Pressure is now in bar, everywhere. */
453 #define factor(d,m) (ir->compress[d][m]*dt/ir->tau_p)
455 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
456 * necessary for triclinic scaling
463 mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
466 case epctSEMIISOTROPIC:
468 mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
470 1.0 - factor(ZZ,ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
472 case epctANISOTROPIC:
475 mu[d][n] = (d==n ? 1.0 : 0.0)
476 -factor(d,n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
478 case epctSURFACETENSION:
479 /* ir->ref_p[0/1] is the reference surface-tension times *
480 * the number of surfaces */
481 if (ir->compress[ZZ][ZZ])
482 p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
484 /* when the compressibity is zero, set the pressure correction *
485 * in the z-direction to zero to get the correct surface tension */
487 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
488 for(d=0; d<DIM-1; d++)
489 mu[d][d] = 1.0 + factor(d,d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
490 - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
493 gmx_fatal(FARGS,"Berendsen pressure coupling type %s not supported yet\n",
494 EPCOUPLTYPETYPE(ir->epct));
497 /* To fullfill the orientation restrictions on triclinic boxes
498 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
499 * the other elements of mu to first order.
501 mu[YY][XX] += mu[XX][YY];
502 mu[ZZ][XX] += mu[XX][ZZ];
503 mu[ZZ][YY] += mu[YY][ZZ];
509 pr_rvecs(debug,0,"PC: pres ",pres,3);
510 pr_rvecs(debug,0,"PC: mu ",mu,3);
513 if (mu[XX][XX]<0.99 || mu[XX][XX]>1.01 ||
514 mu[YY][YY]<0.99 || mu[YY][YY]>1.01 ||
515 mu[ZZ][ZZ]<0.99 || mu[ZZ][ZZ]>1.01) {
517 sprintf(buf,"\nStep %s Warning: pressure scaling more than 1%%, "
519 gmx_step_str(step,buf2),mu[XX][XX],mu[YY][YY],mu[ZZ][ZZ]);
521 fprintf(fplog,"%s",buf);
522 fprintf(stderr,"%s",buf);
526 void berendsen_pscale(t_inputrec *ir,matrix mu,
527 matrix box,matrix box_rel,
528 int start,int nr_atoms,
529 rvec x[],unsigned short cFREEZE[],
532 ivec *nFreeze=ir->opts.nFreeze;
535 /* Scale the positions */
536 for (n=start; n<start+nr_atoms; n++) {
541 x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
543 x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
545 x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
547 /* compute final boxlengths */
548 for (d=0; d<DIM; d++) {
549 box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
550 box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
551 box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
554 preserve_box_shape(ir,box_rel,box);
556 /* (un)shifting should NOT be done after this,
557 * since the box vectors might have changed
559 inc_nrnb(nrnb,eNR_PCOUPL,nr_atoms);
562 void berendsen_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt)
570 for(i=0; (i<opts->ngtc); i++)
574 T = ekind->tcstat[i].T;
578 T = ekind->tcstat[i].Th;
581 if ((opts->tau_t[i] > 0) && (T > 0.0)) {
582 reft = max(0.0,opts->ref_t[i]);
583 lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
584 ekind->tcstat[i].lambda = max(min(lll,1.25),0.8);
587 ekind->tcstat[i].lambda = 1.0;
592 fprintf(debug,"TC: group %d: T: %g, Lambda: %g\n",
593 i,T,ekind->tcstat[i].lambda);
598 static int poisson_variate(real lambda,gmx_rng_t rng) {
609 p *= gmx_rng_uniform_real(rng);
615 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)
618 int i,j,k,d,len,n,ngtc,gc=0;
619 int nshake, nsettle, nrandom, nrand_group;
620 real boltz,scal,reft,prand;
623 /* convenience variables */
627 /* idef is only passed in if it's chance-per-particle andersen, so
628 it essentially serves as a boolean to determine which type of
629 andersen is being used */
632 /* randomly atoms to randomize. However, all constraint
633 groups have to have either all of the atoms or none of the
637 1. Select whether or not to randomize each atom to get the correct probability.
638 2. Cycle through the constraint groups.
639 2a. for each constraint group, determine the fraction f of that constraint group that are
640 chosen to be randomized.
641 2b. all atoms in the constraint group are randomized with probability f.
645 if ((rate < 0.05) && (md->homenr > 50))
647 /* if the rate is relatively high, use a standard method, if low rate,
649 /* poisson distributions approxmation, more efficient for
650 * low rates, fewer random numbers required */
651 nrandom = poisson_variate(md->homenr*rate,rng); /* how many do we randomize? Use Poisson. */
652 /* now we know how many, choose them randomly. No worries about repeats, at this rate, it's negligible.
653 worst thing that happens, it lowers the true rate an negligible amount */
654 for (i=0;i<nrandom;i++)
656 randatom[(int)(gmx_rng_uniform_real(rng)*md->homenr)] = TRUE;
661 for (i=0;i<md->homenr;i++)
663 if (gmx_rng_uniform_real(rng)<rate)
671 /* instead of looping over the constraint groups, if we had a
672 list of which atoms were in which constraint groups, we
673 could then loop over only the groups that are randomized
674 now. But that is not available now. Create later after
675 determining whether there actually is any slowing. */
677 /* first, loop through the settles to make sure all groups either entirely randomized, or not randomized. */
679 nsettle = idef->il[F_SETTLE].nr/2;
680 for (i=0;i<nsettle;i++)
682 iatoms = idef->il[F_SETTLE].iatoms;
684 for (k=0;k<3;k++) /* settles are always 3 atoms, hardcoded */
686 if (randatom[iatoms[2*i+1]+k])
688 nrand_group++; /* count the number of atoms to be shaken in the settles group */
689 randatom[iatoms[2*i+1]+k] = FALSE;
695 prand = (nrand_group)/3.0; /* use this fraction to compute the probability the
696 whole group is randomized */
697 if (gmx_rng_uniform_real(rng)<prand)
701 randatom[iatoms[2*i+1]+k] = TRUE; /* mark them all to be randomized */
708 /* now loop through the shake groups */
710 for (i=0;i<nshake;i++)
712 iatoms = &(idef->il[F_CONSTR].iatoms[sblock[i]]);
713 len = sblock[i+1]-sblock[i];
718 { /* only 2/3 of the sblock items are atoms, the others are labels */
719 if (randatom[iatoms[k]])
722 randatom[iatoms[k]] = FALSE; /* need to mark it false here in case the atom is in more than
723 one group in the shake block */
730 prand = (nrand_group)/(1.0*(2*len/3));
731 if (gmx_rng_uniform_real(rng)<prand)
736 { /* only 2/3 of the sblock items are atoms, the others are labels */
737 randatom[iatoms[k]] = TRUE; /* randomize all of them */
747 for (i=0;i<md->homenr;i++) /* now loop over the list of atoms */
751 randatom_list[n] = i;
755 nrandom = n; /* there are some values of nrandom for which
756 this algorithm won't work; for example all
757 water molecules and nrandom =/= 3. Better to
758 recount and use this number (which we
759 calculate anyway: it will not affect
760 the average number of atoms accepted.
766 /* if it's andersen-massive, then randomize all the atoms */
767 nrandom = md->homenr;
768 for (i=0;i<nrandom;i++)
770 randatom_list[i] = i;
774 /* randomize the velocities of the selected particles */
776 for (i=0;i<nrandom;i++) /* now loop over the list of atoms */
778 n = randatom_list[i];
781 gc = md->cTC[n]; /* assign the atom to a temperature group if there are more than one */
785 scal = sqrt(boltzfac[gc]*md->invmass[n]);
788 state->v[n][d] = scal*gmx_rng_gaussian_table(rng);
791 randatom[n] = FALSE; /* unmark this atom for randomization */
796 void nosehoover_tcoupl(t_grpopts *opts,gmx_ekindata_t *ekind,real dt,
797 double xi[],double vxi[], t_extmass *MassQ)
802 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
804 for(i=0; (i<opts->ngtc); i++)
806 reft = max(0.0,opts->ref_t[i]);
808 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
809 xi[i] += dt*(oldvxi + vxi[i])*0.5;
813 t_state *init_bufstate(const t_state *template_state)
816 int nc = template_state->nhchainlength;
818 snew(state->nosehoover_xi,nc*template_state->ngtc);
819 snew(state->nosehoover_vxi,nc*template_state->ngtc);
820 snew(state->therm_integral,template_state->ngtc);
821 snew(state->nhpres_xi,nc*template_state->nnhpres);
822 snew(state->nhpres_vxi,nc*template_state->nnhpres);
827 void destroy_bufstate(t_state *state)
831 sfree(state->nosehoover_xi);
832 sfree(state->nosehoover_vxi);
833 sfree(state->therm_integral);
834 sfree(state->nhpres_xi);
835 sfree(state->nhpres_vxi);
839 void trotter_update(t_inputrec *ir,gmx_large_int_t step, gmx_ekindata_t *ekind,
840 gmx_enerdata_t *enerd, t_state *state,
841 tensor vir, t_mdatoms *md,
842 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
845 int n,i,j,d,ntgrp,ngtc,gc=0;
846 t_grp_tcstat *tcstat;
848 gmx_large_int_t step_eff;
849 real ecorr,pcorr,dvdlcorr;
850 real bmass,qmass,reft,kT,dt,nd;
851 tensor dumpres,dumvir;
852 double *scalefac,dtc;
854 rvec sumv={0,0,0},consk;
857 if (trotter_seqno <= ettTSEQ2)
859 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
860 is actually the last half step from the previous step. Thus the first half step
861 actually corresponds to the n-1 step*/
867 bCouple = (ir->nsttcouple == 1 ||
868 do_per_step(step_eff+ir->nsttcouple,ir->nsttcouple));
870 trotter_seq = trotter_seqlist[trotter_seqno];
872 /* signal we are returning if nothing is going to be done in this routine */
873 if ((trotter_seq[0] == etrtSKIPALL) || !(bCouple))
878 dtc = ir->nsttcouple*ir->delta_t;
879 opts = &(ir->opts); /* just for ease of referencing */
882 snew(scalefac,opts->ngtc);
887 /* execute the series of trotter updates specified in the trotterpart array */
889 for (i=0;i<NTROTTERPARTS;i++){
890 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
891 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
900 switch (trotter_seq[i])
904 boxv_trotter(ir,&(state->veta),dt,state->box,ekind,vir,
905 enerd->term[F_PDISPCORR],enerd->term[F_DISPCORR],MassQ);
909 NHC_trotter(opts,state->nnhpres,ekind,dt,state->nhpres_xi,
910 state->nhpres_vxi,NULL,&(state->veta),MassQ,FALSE);
914 NHC_trotter(opts,opts->ngtc,ekind,dt,state->nosehoover_xi,
915 state->nosehoover_vxi,scalefac,NULL,MassQ,(ir->eI==eiVV));
916 /* need to rescale the kinetic energies and velocities here. Could
917 scale the velocities later, but we need them scaled in order to
918 produce the correct outputs, so we'll scale them here. */
920 for (i=0; i<ngtc;i++)
922 tcstat = &ekind->tcstat[i];
923 tcstat->vscale_nhc = scalefac[i];
924 tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
925 tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
927 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
928 /* but do we actually need the total? */
930 /* modify the velocities as well */
931 for (n=md->start;n<md->start+md->homenr;n++)
933 if (md->cTC) /* does this conditional need to be here? is this always true?*/
939 state->v[n][d] *= scalefac[gc];
946 sumv[d] += (state->v[n][d])/md->invmass[n];
955 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
963 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
965 fprintf(debug,"Conserved kappa: %15.8f %15.8f %15.8f\n",consk[0],consk[1],consk[2]);
973 extern void init_npt_masses(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bInit)
975 int n,i,j,d,ntgrp,ngtc,nnhpres,nh,gc=0;
976 t_grp_tcstat *tcstat;
978 real ecorr,pcorr,dvdlcorr;
979 real bmass,qmass,reft,kT,dt,ndj,nd;
980 tensor dumpres,dumvir;
982 opts = &(ir->opts); /* just for ease of referencing */
983 ngtc = ir->opts.ngtc;
984 nnhpres = state->nnhpres;
985 nh = state->nhchainlength;
987 if (ir->eI == eiMD) {
990 snew(MassQ->Qinv,ngtc);
992 for(i=0; (i<ngtc); i++)
994 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
996 MassQ->Qinv[i]=1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
1004 else if (EI_VV(ir->eI))
1006 /* Set pressure variables */
1010 if (state->vol0 == 0)
1012 state->vol0 = det(state->box);
1013 /* because we start by defining a fixed
1014 compressibility, we need the volume at this
1015 compressibility to solve the problem. */
1019 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
1020 /* Consider evaluating eventually if this the right mass to use. All are correct, some might be more stable */
1021 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
1022 /* An alternate mass definition, from Tuckerman et al. */
1023 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
1028 MassQ->Winvm[d][n]= PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
1029 /* not clear this is correct yet for the anisotropic case. Will need to reevaluate
1030 before using MTTK for anisotropic states.*/
1033 /* Allocate space for thermostat variables */
1036 snew(MassQ->Qinv,ngtc*nh);
1039 /* now, set temperature variables */
1040 for (i=0; i<ngtc; i++)
1042 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
1044 reft = max(0.0,opts->ref_t[i]);
1057 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
1065 MassQ->Qinv[i*nh+j] = 0.0;
1072 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
1074 int n,i,j,d,ntgrp,ngtc,nnhpres,nh,gc=0;
1075 t_grp_tcstat *tcstat;
1077 real ecorr,pcorr,dvdlcorr;
1078 real bmass,qmass,reft,kT,dt,ndj,nd;
1079 tensor dumpres,dumvir;
1082 opts = &(ir->opts); /* just for ease of referencing */
1084 nnhpres = state->nnhpres;
1085 nh = state->nhchainlength;
1087 init_npt_masses(ir, state, MassQ, TRUE);
1089 /* first, initialize clear all the trotter calls */
1090 snew(trotter_seq,ettTSEQMAX);
1091 for (i=0;i<ettTSEQMAX;i++)
1093 snew(trotter_seq[i],NTROTTERPARTS);
1094 for (j=0;j<NTROTTERPARTS;j++) {
1095 trotter_seq[i][j] = etrtNONE;
1097 trotter_seq[i][0] = etrtSKIPALL;
1102 /* no trotter calls, so we never use the values in the array.
1103 * We access them (so we need to define them, but ignore
1109 /* compute the kinetic energy by using the half step velocities or
1110 * the kinetic energies, depending on the order of the trotter calls */
1114 if (IR_NPT_TROTTER(ir))
1116 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1117 We start with the initial one. */
1118 /* first, a round that estimates veta. */
1119 trotter_seq[0][0] = etrtBAROV;
1121 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1123 /* The first half trotter update */
1124 trotter_seq[2][0] = etrtBAROV;
1125 trotter_seq[2][1] = etrtNHC;
1126 trotter_seq[2][2] = etrtBARONHC;
1128 /* The second half trotter update */
1129 trotter_seq[3][0] = etrtBARONHC;
1130 trotter_seq[3][1] = etrtNHC;
1131 trotter_seq[3][2] = etrtBAROV;
1133 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1136 else if (IR_NVT_TROTTER(ir))
1138 /* This is the easy version - there are only two calls, both the same.
1139 Otherwise, even easier -- no calls */
1140 trotter_seq[2][0] = etrtNHC;
1141 trotter_seq[3][0] = etrtNHC;
1143 else if (IR_NPH_TROTTER(ir))
1145 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1146 We start with the initial one. */
1147 /* first, a round that estimates veta. */
1148 trotter_seq[0][0] = etrtBAROV;
1150 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
1152 /* The first half trotter update */
1153 trotter_seq[2][0] = etrtBAROV;
1154 trotter_seq[2][1] = etrtBARONHC;
1156 /* The second half trotter update */
1157 trotter_seq[3][0] = etrtBARONHC;
1158 trotter_seq[3][1] = etrtBAROV;
1160 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
1163 else if (ir->eI==eiVVAK)
1165 if (IR_NPT_TROTTER(ir))
1167 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1168 We start with the initial one. */
1169 /* first, a round that estimates veta. */
1170 trotter_seq[0][0] = etrtBAROV;
1172 /* The first half trotter update, part 1 -- double update, because it commutes */
1173 trotter_seq[1][0] = etrtNHC;
1175 /* The first half trotter update, part 2 */
1176 trotter_seq[2][0] = etrtBAROV;
1177 trotter_seq[2][1] = etrtBARONHC;
1179 /* The second half trotter update, part 1 */
1180 trotter_seq[3][0] = etrtBARONHC;
1181 trotter_seq[3][1] = etrtBAROV;
1183 /* The second half trotter update */
1184 trotter_seq[4][0] = etrtNHC;
1186 else if (IR_NVT_TROTTER(ir))
1188 /* This is the easy version - there is only one call, both the same.
1189 Otherwise, even easier -- no calls */
1190 trotter_seq[1][0] = etrtNHC;
1191 trotter_seq[4][0] = etrtNHC;
1193 else if (IR_NPH_TROTTER(ir))
1195 /* This is the complicated version - there are 4 possible calls, depending on ordering.
1196 We start with the initial one. */
1197 /* first, a round that estimates veta. */
1198 trotter_seq[0][0] = etrtBAROV;
1200 /* The first half trotter update, part 1 -- leave zero */
1201 trotter_seq[1][0] = etrtNHC;
1203 /* The first half trotter update, part 2 */
1204 trotter_seq[2][0] = etrtBAROV;
1205 trotter_seq[2][1] = etrtBARONHC;
1207 /* The second half trotter update, part 1 */
1208 trotter_seq[3][0] = etrtBARONHC;
1209 trotter_seq[3][1] = etrtBAROV;
1211 /* The second half trotter update -- blank for now */
1219 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
1222 snew(MassQ->QPinv,nnhpres*opts->nhchainlength);
1224 /* barostat temperature */
1225 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1227 reft = max(0.0,opts->ref_t[0]);
1229 for (i=0;i<nnhpres;i++) {
1239 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1245 for (i=0;i<nnhpres;i++) {
1248 MassQ->QPinv[i*nh+j]=0.0;
1255 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1257 int i,j,nd,ndj,bmass,qmass,ngtcall;
1258 real ener_npt,reft,eta,kT,tau;
1262 int nh = state->nhchainlength;
1266 /* now we compute the contribution of the pressure to the conserved quantity*/
1268 if (ir->epc==epcMTTK)
1270 /* find the volume, and the kinetic energy of the volume */
1275 /* contribution from the pressure momenenta */
1276 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1278 /* contribution from the PV term */
1279 vol = det(state->box);
1280 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1283 case epctANISOTROPIC:
1287 case epctSURFACETENSION:
1290 case epctSEMIISOTROPIC:
1298 if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1300 /* add the energy from the barostat thermostat chain */
1301 for (i=0;i<state->nnhpres;i++) {
1303 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1304 ivxi = &state->nhpres_vxi[i*nh];
1305 ixi = &state->nhpres_xi[i*nh];
1306 iQinv = &(MassQ->QPinv[i*nh]);
1307 reft = max(ir->opts.ref_t[0],0); /* using 'System' temperature */
1314 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1315 /* contribution from the thermal variable of the NH chain */
1316 ener_npt += ixi[j]*kT;
1320 fprintf(debug,"P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f",i,j,ivxi[j],ixi[j]);
1328 for(i=0; i<ir->opts.ngtc; i++)
1330 ixi = &state->nosehoover_xi[i*nh];
1331 ivxi = &state->nosehoover_vxi[i*nh];
1332 iQinv = &(MassQ->Qinv[i*nh]);
1334 nd = ir->opts.nrdf[i];
1335 reft = max(ir->opts.ref_t[i],0);
1340 if (IR_NVT_TROTTER(ir))
1342 /* contribution from the thermal momenta of the NH chain */
1346 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1347 /* contribution from the thermal variable of the NH chain */
1355 ener_npt += ndj*ixi[j]*kT;
1359 else /* Other non Trotter temperature NH control -- no chains yet. */
1361 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1362 ener_npt += nd*ixi[0]*kT;
1370 static real vrescale_gamdev(int ia, gmx_rng_t rng)
1371 /* Gamma distribution, adapted from numerical recipes */
1374 real am,e,s,v1,v2,x,y;
1381 for(j=1; j<=ia; j++)
1383 x *= gmx_rng_uniform_real(rng);
1397 v1 = gmx_rng_uniform_real(rng);
1398 v2 = 2.0*gmx_rng_uniform_real(rng)-1.0;
1400 while (v1*v1 + v2*v2 > 1.0 ||
1401 v1*v1*GMX_REAL_MAX < 3.0*ia);
1402 /* The last check above ensures that both x (3.0 > 2.0 in s)
1403 * and the pre-factor for e do not go out of range.
1407 s = sqrt(2.0*am + 1.0);
1411 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1413 while (gmx_rng_uniform_real(rng) > e);
1419 static real vrescale_sumnoises(int nn,gmx_rng_t rng)
1422 * Returns the sum of n independent gaussian noises squared
1423 * (i.e. equivalent to summing the square of the return values
1424 * of nn calls to gmx_rng_gaussian_real).xs
1430 } else if (nn == 1) {
1431 rr = gmx_rng_gaussian_real(rng);
1433 } else if (nn % 2 == 0) {
1434 return 2.0*vrescale_gamdev(nn/2,rng);
1436 rr = gmx_rng_gaussian_real(rng);
1437 return 2.0*vrescale_gamdev((nn-1)/2,rng) + rr*rr;
1441 static real vrescale_resamplekin(real kk,real sigma, int ndeg, real taut,
1445 * Generates a new value for the kinetic energy,
1446 * according to Bussi et al JCP (2007), Eq. (A7)
1447 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1448 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1449 * ndeg: number of degrees of freedom of the atoms to be thermalized
1450 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1455 factor = exp(-1.0/taut);
1459 rr = gmx_rng_gaussian_real(rng);
1462 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1,rng) + rr*rr)/ndeg - kk) +
1463 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1466 void vrescale_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt,
1467 double therm_integral[],gmx_rng_t rng)
1471 real Ek,Ek_ref1,Ek_ref,Ek_new;
1475 for(i=0; (i<opts->ngtc); i++)
1479 Ek = trace(ekind->tcstat[i].ekinf);
1483 Ek = trace(ekind->tcstat[i].ekinh);
1486 if (opts->tau_t[i] > 0 && opts->nrdf[i] > 0 && Ek > 0)
1488 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1489 Ek_ref = Ek_ref1*opts->nrdf[i];
1491 Ek_new = vrescale_resamplekin(Ek,Ek_ref,opts->nrdf[i],
1492 opts->tau_t[i]/dt,rng);
1494 /* Analytically Ek_new>=0, but we check for rounding errors */
1497 ekind->tcstat[i].lambda = 0.0;
1501 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1504 therm_integral[i] -= Ek_new - Ek;
1508 fprintf(debug,"TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1509 i,Ek_ref,Ek,Ek_new,ekind->tcstat[i].lambda);
1514 ekind->tcstat[i].lambda = 1.0;
1519 real vrescale_energy(t_grpopts *opts,double therm_integral[])
1525 for(i=0; i<opts->ngtc; i++) {
1526 ener += therm_integral[i];
1532 void rescale_velocities(gmx_ekindata_t *ekind,t_mdatoms *mdatoms,
1533 int start,int end,rvec v[])
1536 t_grp_tcstat *tcstat;
1537 unsigned short *cACC,*cTC;
1542 tcstat = ekind->tcstat;
1547 gstat = ekind->grpstat;
1548 cACC = mdatoms->cACC;
1552 for(n=start; n<end; n++)
1562 /* Only scale the velocity component relative to the COM velocity */
1563 rvec_sub(v[n],gstat[ga].u,vrel);
1564 lg = tcstat[gt].lambda;
1565 for(d=0; d<DIM; d++)
1567 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1574 for(n=start; n<end; n++)
1580 lg = tcstat[gt].lambda;
1581 for(d=0; d<DIM; d++)
1590 /* set target temperatures if we are annealing */
1591 void update_annealing_target_temp(t_grpopts *opts,real t)
1594 real pert,thist=0,x;
1596 for(i=0;i<opts->ngtc;i++) {
1597 npoints = opts->anneal_npoints[i];
1598 switch (opts->annealing[i]) {
1602 /* calculate time modulo the period */
1603 pert = opts->anneal_time[i][npoints-1];
1605 thist = t - n*pert; /* modulo time */
1606 /* Make sure rounding didn't get us outside the interval */
1607 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1614 gmx_fatal(FARGS,"Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",i,opts->ngtc,npoints);
1616 /* We are doing annealing for this group if we got here,
1617 * and we have the (relative) time as thist.
1618 * calculate target temp */
1620 while ((j < npoints-1) && (thist>(opts->anneal_time[i][j+1])))
1622 if (j < npoints-1) {
1623 /* Found our position between points j and j+1.
1624 * Interpolate: x is the amount from j+1, (1-x) from point j
1625 * First treat possible jumps in temperature as a special case.
1627 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1628 opts->ref_t[i]=opts->anneal_temp[i][j+1];
1630 x = ((thist-opts->anneal_time[i][j])/
1631 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1632 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1636 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];