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; /* 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);
247 pscal = calc_pres(ir->ePBC,nwall,box,ekinmod,vir,localpres) + pcorr;
250 GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
256 * This file implements temperature and pressure coupling algorithms:
257 * For now only the Weak coupling and the modified weak coupling.
259 * Furthermore computation of pressure and temperature is done here
263 real calc_pres(int ePBC,int nwall,matrix box,tensor ekin,tensor vir,
269 if (ePBC==epbcNONE || (ePBC==epbcXY && nwall!=2))
272 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
273 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
277 fac=PRESFAC*2.0/det(box);
278 for(n=0; (n<DIM); n++)
279 for(m=0; (m<DIM); m++)
280 pres[n][m] = (ekin[n][m] - vir[n][m])*fac;
283 pr_rvecs(debug,0,"PC: pres",pres,DIM);
284 pr_rvecs(debug,0,"PC: ekin",ekin,DIM);
285 pr_rvecs(debug,0,"PC: vir ",vir, DIM);
286 pr_rvecs(debug,0,"PC: box ",box, DIM);
289 return trace(pres)/DIM;
292 real calc_temp(real ekin,real nrdf)
295 return (2.0*ekin)/(nrdf*BOLTZ);
300 void parrinellorahman_pcoupl(FILE *fplog,gmx_large_int_t step,
301 t_inputrec *ir,real dt,tensor pres,
302 tensor box,tensor box_rel,tensor boxv,
303 tensor M,matrix mu,gmx_bool bFirstStep)
305 /* This doesn't do any coordinate updating. It just
306 * integrates the box vector equations from the calculated
307 * acceleration due to pressure difference. We also compute
308 * the tensor M which is used in update to couple the particle
309 * coordinates to the box vectors.
311 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
314 * M_nk = (h') * (h' * h + h' h) * h
316 * with the dots denoting time derivatives and h is the transformation from
317 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
318 * This also goes for the pressure and M tensors - they are transposed relative
319 * to ours. Our equation thus becomes:
322 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
324 * where b is the gromacs box matrix.
325 * Our box accelerations are given by
327 * b = vol/W inv(box') * (P-ref_P) (=h')
332 real vol=box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
333 real atot,arel,change,maxchange,xy_pressure;
334 tensor invbox,pdiff,t1,t2;
338 m_inv_ur0(box,invbox);
341 /* Note that PRESFAC does not occur here.
342 * The pressure and compressibility always occur as a product,
343 * therefore the pressure unit drops out.
345 maxl=max(box[XX][XX],box[YY][YY]);
346 maxl=max(maxl,box[ZZ][ZZ]);
350 (4*M_PI*M_PI*ir->compress[d][n])/(3*ir->tau_p*ir->tau_p*maxl);
352 m_sub(pres,ir->ref_p,pdiff);
354 if(ir->epct==epctSURFACETENSION) {
355 /* Unlike Berendsen coupling it might not be trivial to include a z
356 * pressure correction here? On the other hand we don't scale the
357 * box momentarily, but change accelerations, so it might not be crucial.
359 xy_pressure=0.5*(pres[XX][XX]+pres[YY][YY]);
361 pdiff[d][d]=(xy_pressure-(pres[ZZ][ZZ]-ir->ref_p[d][d]/box[d][d]));
364 tmmul(invbox,pdiff,t1);
365 /* Move the off-diagonal elements of the 'force' to one side to ensure
366 * that we obey the box constraints.
370 t1[d][n] += t1[n][d];
376 case epctANISOTROPIC:
379 t1[d][n] *= winv[d][n]*vol;
382 /* calculate total volume acceleration */
383 atot=box[XX][XX]*box[YY][YY]*t1[ZZ][ZZ]+
384 box[XX][XX]*t1[YY][YY]*box[ZZ][ZZ]+
385 t1[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
387 /* set all RELATIVE box accelerations equal, and maintain total V
391 t1[d][n] = winv[0][0]*vol*arel*box[d][n];
393 case epctSEMIISOTROPIC:
394 case epctSURFACETENSION:
395 /* Note the correction to pdiff above for surftens. coupling */
397 /* calculate total XY volume acceleration */
398 atot=box[XX][XX]*t1[YY][YY]+t1[XX][XX]*box[YY][YY];
399 arel=atot/(2*box[XX][XX]*box[YY][YY]);
400 /* set RELATIVE XY box accelerations equal, and maintain total V
401 * change speed. Dont change the third box vector accelerations */
404 t1[d][n] = winv[d][n]*vol*arel*box[d][n];
406 t1[ZZ][n] *= winv[d][n]*vol;
409 gmx_fatal(FARGS,"Parrinello-Rahman pressure coupling type %s "
410 "not supported yet\n",EPCOUPLTYPETYPE(ir->epct));
417 boxv[d][n] += dt*t1[d][n];
419 /* We do NOT update the box vectors themselves here, since
420 * we need them for shifting later. It is instead done last
421 * in the update() routine.
424 /* Calculate the change relative to diagonal elements-
425 since it's perfectly ok for the off-diagonal ones to
426 be zero it doesn't make sense to check the change relative
430 change=fabs(dt*boxv[d][n]/box[d][d]);
432 if (change>maxchange)
436 if (maxchange > 0.01 && fplog) {
439 "\nStep %s Warning: Pressure scaling more than 1%%. "
440 "This may mean your system\n is not yet equilibrated. "
441 "Use of Parrinello-Rahman pressure coupling during\n"
442 "equilibration can lead to simulation instability, "
443 "and is discouraged.\n",
444 gmx_step_str(step,buf));
448 preserve_box_shape(ir,box_rel,boxv);
450 mtmul(boxv,box,t1); /* t1=boxv * b' */
454 /* Determine the scaling matrix mu for the coordinates */
457 t1[d][n] = box[d][n] + dt*boxv[d][n];
458 preserve_box_shape(ir,box_rel,t1);
459 /* t1 is the box at t+dt, determine mu as the relative change */
460 mmul_ur0(invbox,t1,mu);
463 void berendsen_pcoupl(FILE *fplog,gmx_large_int_t step,
464 t_inputrec *ir,real dt, tensor pres,matrix box,
468 real scalar_pressure, xy_pressure, p_corr_z;
469 char *ptr,buf[STRLEN];
472 * Calculate the scaling matrix mu
476 for(d=0; d<DIM; d++) {
477 scalar_pressure += pres[d][d]/DIM;
479 xy_pressure += pres[d][d]/(DIM-1);
481 /* Pressure is now in bar, everywhere. */
482 #define factor(d,m) (ir->compress[d][m]*dt/ir->tau_p)
484 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
485 * necessary for triclinic scaling
492 mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
495 case epctSEMIISOTROPIC:
497 mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
499 1.0 - factor(ZZ,ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
501 case epctANISOTROPIC:
504 mu[d][n] = (d==n ? 1.0 : 0.0)
505 -factor(d,n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
507 case epctSURFACETENSION:
508 /* ir->ref_p[0/1] is the reference surface-tension times *
509 * the number of surfaces */
510 if (ir->compress[ZZ][ZZ])
511 p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
513 /* when the compressibity is zero, set the pressure correction *
514 * in the z-direction to zero to get the correct surface tension */
516 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
517 for(d=0; d<DIM-1; d++)
518 mu[d][d] = 1.0 + factor(d,d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
519 - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
522 gmx_fatal(FARGS,"Berendsen pressure coupling type %s not supported yet\n",
523 EPCOUPLTYPETYPE(ir->epct));
526 /* To fullfill the orientation restrictions on triclinic boxes
527 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
528 * the other elements of mu to first order.
530 mu[YY][XX] += mu[XX][YY];
531 mu[ZZ][XX] += mu[XX][ZZ];
532 mu[ZZ][YY] += mu[YY][ZZ];
538 pr_rvecs(debug,0,"PC: pres ",pres,3);
539 pr_rvecs(debug,0,"PC: mu ",mu,3);
542 if (mu[XX][XX]<0.99 || mu[XX][XX]>1.01 ||
543 mu[YY][YY]<0.99 || mu[YY][YY]>1.01 ||
544 mu[ZZ][ZZ]<0.99 || mu[ZZ][ZZ]>1.01) {
546 sprintf(buf,"\nStep %s Warning: pressure scaling more than 1%%, "
548 gmx_step_str(step,buf2),mu[XX][XX],mu[YY][YY],mu[ZZ][ZZ]);
550 fprintf(fplog,"%s",buf);
551 fprintf(stderr,"%s",buf);
555 void berendsen_pscale(t_inputrec *ir,matrix mu,
556 matrix box,matrix box_rel,
557 int start,int nr_atoms,
558 rvec x[],unsigned short cFREEZE[],
561 ivec *nFreeze=ir->opts.nFreeze;
564 /* Scale the positions */
565 for (n=start; n<start+nr_atoms; n++) {
570 x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
572 x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
574 x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
576 /* compute final boxlengths */
577 for (d=0; d<DIM; d++) {
578 box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
579 box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
580 box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
583 preserve_box_shape(ir,box_rel,box);
585 /* (un)shifting should NOT be done after this,
586 * since the box vectors might have changed
588 inc_nrnb(nrnb,eNR_PCOUPL,nr_atoms);
591 void berendsen_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt)
599 for(i=0; (i<opts->ngtc); i++)
603 T = ekind->tcstat[i].T;
607 T = ekind->tcstat[i].Th;
610 if ((opts->tau_t[i] > 0) && (T > 0.0)) {
612 reft = max(0.0,opts->ref_t[i]);
613 lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
614 ekind->tcstat[i].lambda = max(min(lll,1.25),0.8);
617 ekind->tcstat[i].lambda = 1.0;
621 fprintf(debug,"TC: group %d: T: %g, Lambda: %g\n",
622 i,T,ekind->tcstat[i].lambda);
626 void nosehoover_tcoupl(t_grpopts *opts,gmx_ekindata_t *ekind,real dt,
627 double xi[],double vxi[], t_extmass *MassQ)
632 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
634 for(i=0; (i<opts->ngtc); i++) {
635 reft = max(0.0,opts->ref_t[i]);
637 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
638 xi[i] += dt*(oldvxi + vxi[i])*0.5;
642 t_state *init_bufstate(const t_state *template_state)
645 int nc = template_state->nhchainlength;
647 snew(state->nosehoover_xi,nc*template_state->ngtc);
648 snew(state->nosehoover_vxi,nc*template_state->ngtc);
649 snew(state->therm_integral,template_state->ngtc);
650 snew(state->nhpres_xi,nc*template_state->nnhpres);
651 snew(state->nhpres_vxi,nc*template_state->nnhpres);
656 void destroy_bufstate(t_state *state)
660 sfree(state->nosehoover_xi);
661 sfree(state->nosehoover_vxi);
662 sfree(state->therm_integral);
663 sfree(state->nhpres_xi);
664 sfree(state->nhpres_vxi);
668 void trotter_update(t_inputrec *ir,gmx_large_int_t step, gmx_ekindata_t *ekind,
669 gmx_enerdata_t *enerd, t_state *state,
670 tensor vir, t_mdatoms *md,
671 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
674 int n,i,j,d,ntgrp,ngtc,gc=0;
675 t_grp_tcstat *tcstat;
677 gmx_large_int_t step_eff;
678 real ecorr,pcorr,dvdlcorr;
679 real bmass,qmass,reft,kT,dt,nd;
680 tensor dumpres,dumvir;
681 double *scalefac,dtc;
683 rvec sumv={0,0,0},consk;
686 if (trotter_seqno <= ettTSEQ2)
688 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
689 is actually the last half step from the previous step. Thus the first half step
690 actually corresponds to the n-1 step*/
696 bCouple = (ir->nsttcouple == 1 ||
697 do_per_step(step_eff+ir->nsttcouple,ir->nsttcouple));
699 trotter_seq = trotter_seqlist[trotter_seqno];
701 /* signal we are returning if nothing is going to be done in this routine */
702 if ((trotter_seq[0] == etrtSKIPALL) || !(bCouple))
707 dtc = ir->nsttcouple*ir->delta_t;
708 opts = &(ir->opts); /* just for ease of referencing */
711 snew(scalefac,opts->ngtc);
716 /* execute the series of trotter updates specified in the trotterpart array */
718 for (i=0;i<NTROTTERPARTS;i++){
719 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
720 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
729 switch (trotter_seq[i])
733 boxv_trotter(ir,&(state->veta),dt,state->box,ekind,vir,
734 enerd->term[F_PDISPCORR],enerd->term[F_DISPCORR],MassQ);
738 NHC_trotter(opts,state->nnhpres,ekind,dt,state->nhpres_xi,
739 state->nhpres_vxi,NULL,&(state->veta),MassQ,FALSE);
743 NHC_trotter(opts,opts->ngtc,ekind,dt,state->nosehoover_xi,
744 state->nosehoover_vxi,scalefac,NULL,MassQ,(ir->eI==eiVV));
745 /* need to rescale the kinetic energies and velocities here. Could
746 scale the velocities later, but we need them scaled in order to
747 produce the correct outputs, so we'll scale them here. */
749 for (i=0; i<ngtc;i++)
751 tcstat = &ekind->tcstat[i];
752 tcstat->vscale_nhc = scalefac[i];
753 tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
754 tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
756 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
757 /* but do we actually need the total? */
759 /* modify the velocities as well */
760 for (n=md->start;n<md->start+md->homenr;n++)
768 state->v[n][d] *= scalefac[gc];
775 sumv[d] += (state->v[n][d])/md->invmass[n];
784 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
792 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
794 fprintf(debug,"Conserved kappa: %15.8f %15.8f %15.8f\n",consk[0],consk[1],consk[2]);
801 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
803 int n,i,j,d,ntgrp,ngtc,nnhpres,nh,gc=0;
804 t_grp_tcstat *tcstat;
806 real ecorr,pcorr,dvdlcorr;
807 real bmass,qmass,reft,kT,dt,ndj,nd;
808 tensor dumpres,dumvir;
811 opts = &(ir->opts); /* just for ease of referencing */
813 nnhpres = state->nnhpres;
814 nh = state->nhchainlength;
816 if (ir->eI == eiMD) {
817 snew(MassQ->Qinv,ngtc);
818 for(i=0; (i<ngtc); i++)
820 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
822 MassQ->Qinv[i]=1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
830 else if (EI_VV(ir->eI))
832 /* Set pressure variables */
834 if (state->vol0 == 0)
836 state->vol0 = det(state->box); /* because we start by defining a fixed compressibility,
837 we need the volume at this compressibility to solve the problem */
840 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
841 /* Investigate this more -- is this the right mass to make it? */
842 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
843 /* An alternate mass definition, from Tuckerman et al. */
844 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
849 MassQ->Winvm[d][n]= PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
850 /* not clear this is correct yet for the anisotropic case*/
853 /* Allocate space for thermostat variables */
854 snew(MassQ->Qinv,ngtc*nh);
856 /* now, set temperature variables */
857 for(i=0; i<ngtc; i++)
859 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
861 reft = max(0.0,opts->ref_t[i]);
874 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
882 MassQ->Qinv[i*nh+j] = 0.0;
888 /* first, initialize clear all the trotter calls */
889 snew(trotter_seq,ettTSEQMAX);
890 for (i=0;i<ettTSEQMAX;i++)
892 snew(trotter_seq[i],NTROTTERPARTS);
893 for (j=0;j<NTROTTERPARTS;j++) {
894 trotter_seq[i][j] = etrtNONE;
896 trotter_seq[i][0] = etrtSKIPALL;
901 /* no trotter calls, so we never use the values in the array.
902 * We access them (so we need to define them, but ignore
908 /* compute the kinetic energy by using the half step velocities or
909 * the kinetic energies, depending on the order of the trotter calls */
913 if (IR_NPT_TROTTER(ir))
915 /* This is the complicated version - there are 4 possible calls, depending on ordering.
916 We start with the initial one. */
917 /* first, a round that estimates veta. */
918 trotter_seq[0][0] = etrtBAROV;
920 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
922 /* The first half trotter update */
923 trotter_seq[2][0] = etrtBAROV;
924 trotter_seq[2][1] = etrtNHC;
925 trotter_seq[2][2] = etrtBARONHC;
927 /* The second half trotter update */
928 trotter_seq[3][0] = etrtBARONHC;
929 trotter_seq[3][1] = etrtNHC;
930 trotter_seq[3][2] = etrtBAROV;
932 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
937 if (IR_NVT_TROTTER(ir))
939 /* This is the easy version - there are only two calls, both the same.
940 Otherwise, even easier -- no calls */
941 trotter_seq[2][0] = etrtNHC;
942 trotter_seq[3][0] = etrtNHC;
945 } else if (ir->eI==eiVVAK) {
946 if (IR_NPT_TROTTER(ir))
948 /* This is the complicated version - there are 4 possible calls, depending on ordering.
949 We start with the initial one. */
950 /* first, a round that estimates veta. */
951 trotter_seq[0][0] = etrtBAROV;
953 /* The first half trotter update, part 1 -- double update, because it commutes */
954 trotter_seq[1][0] = etrtNHC;
956 /* The first half trotter update, part 2 */
957 trotter_seq[2][0] = etrtBAROV;
958 trotter_seq[2][1] = etrtBARONHC;
960 /* The second half trotter update, part 1 */
961 trotter_seq[3][0] = etrtBARONHC;
962 trotter_seq[3][1] = etrtBAROV;
964 /* The second half trotter update -- blank for now */
965 trotter_seq[4][0] = etrtNHC;
969 if (IR_NVT_TROTTER(ir))
971 /* This is the easy version - there is only one call, both the same.
972 Otherwise, even easier -- no calls */
973 trotter_seq[1][0] = etrtNHC;
974 trotter_seq[4][0] = etrtNHC;
983 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
986 snew(MassQ->QPinv,nnhpres*opts->nhchainlength);
988 /* barostat temperature */
989 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
991 reft = max(0.0,opts->ref_t[0]);
993 for (i=0;i<nnhpres;i++) {
1003 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1009 for (i=0;i<nnhpres;i++) {
1012 MassQ->QPinv[i*nh+j]=0.0;
1019 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1021 int i,j,nd,ndj,bmass,qmass,ngtcall;
1022 real ener_npt,reft,eta,kT,tau;
1026 int nh = state->nhchainlength;
1030 /* now we compute the contribution of the pressure to the conserved quantity*/
1032 if (ir->epc==epcMTTK)
1034 /* find the volume, and the kinetic energy of the volume */
1039 /* contribution from the pressure momenenta */
1040 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1042 /* contribution from the PV term */
1043 vol = det(state->box);
1044 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1047 case epctANISOTROPIC:
1051 case epctSURFACETENSION:
1054 case epctSEMIISOTROPIC:
1062 if (IR_NPT_TROTTER(ir))
1064 /* add the energy from the barostat thermostat chain */
1065 for (i=0;i<state->nnhpres;i++) {
1067 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1068 ivxi = &state->nhpres_vxi[i*nh];
1069 ixi = &state->nhpres_xi[i*nh];
1070 iQinv = &(MassQ->QPinv[i*nh]);
1071 reft = max(ir->opts.ref_t[0],0); /* using 'System' temperature */
1078 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1079 /* contribution from the thermal variable of the NH chain */
1080 ener_npt += ixi[j]*kT;
1084 fprintf(debug,"P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f",i,j,ivxi[j],ixi[j]);
1092 for(i=0; i<ir->opts.ngtc; i++)
1094 ixi = &state->nosehoover_xi[i*nh];
1095 ivxi = &state->nosehoover_vxi[i*nh];
1096 iQinv = &(MassQ->Qinv[i*nh]);
1098 nd = ir->opts.nrdf[i];
1099 reft = max(ir->opts.ref_t[i],0);
1104 if (IR_NVT_TROTTER(ir))
1106 /* contribution from the thermal momenta of the NH chain */
1110 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1111 /* contribution from the thermal variable of the NH chain */
1119 ener_npt += ndj*ixi[j]*kT;
1123 else /* Other non Trotter temperature NH control -- no chains yet. */
1125 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1126 ener_npt += nd*ixi[0]*kT;
1134 static real vrescale_gamdev(int ia, gmx_rng_t rng)
1135 /* Gamma distribution, adapted from numerical recipes */
1138 real am,e,s,v1,v2,x,y;
1145 for(j=1; j<=ia; j++)
1147 x *= gmx_rng_uniform_real(rng);
1161 v1 = gmx_rng_uniform_real(rng);
1162 v2 = 2.0*gmx_rng_uniform_real(rng)-1.0;
1164 while (v1*v1 + v2*v2 > 1.0 ||
1165 v1*v1*GMX_REAL_MAX < 3.0*ia);
1166 /* The last check above ensures that both x (3.0 > 2.0 in s)
1167 * and the pre-factor for e do not go out of range.
1171 s = sqrt(2.0*am + 1.0);
1175 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1177 while (gmx_rng_uniform_real(rng) > e);
1183 static real vrescale_sumnoises(int nn,gmx_rng_t rng)
1186 * Returns the sum of n independent gaussian noises squared
1187 * (i.e. equivalent to summing the square of the return values
1188 * of nn calls to gmx_rng_gaussian_real).xs
1194 } else if (nn == 1) {
1195 rr = gmx_rng_gaussian_real(rng);
1197 } else if (nn % 2 == 0) {
1198 return 2.0*vrescale_gamdev(nn/2,rng);
1200 rr = gmx_rng_gaussian_real(rng);
1201 return 2.0*vrescale_gamdev((nn-1)/2,rng) + rr*rr;
1205 static real vrescale_resamplekin(real kk,real sigma, int ndeg, real taut,
1209 * Generates a new value for the kinetic energy,
1210 * according to Bussi et al JCP (2007), Eq. (A7)
1211 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1212 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1213 * ndeg: number of degrees of freedom of the atoms to be thermalized
1214 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1219 factor = exp(-1.0/taut);
1223 rr = gmx_rng_gaussian_real(rng);
1226 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1,rng) + rr*rr)/ndeg - kk) +
1227 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1230 void vrescale_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt,
1231 double therm_integral[],gmx_rng_t rng)
1235 real Ek,Ek_ref1,Ek_ref,Ek_new;
1239 for(i=0; (i<opts->ngtc); i++)
1243 Ek = trace(ekind->tcstat[i].ekinf);
1247 Ek = trace(ekind->tcstat[i].ekinh);
1250 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1252 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1253 Ek_ref = Ek_ref1*opts->nrdf[i];
1255 Ek_new = vrescale_resamplekin(Ek,Ek_ref,opts->nrdf[i],
1256 opts->tau_t[i]/dt,rng);
1258 /* Analytically Ek_new>=0, but we check for rounding errors */
1261 ekind->tcstat[i].lambda = 0.0;
1265 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1268 therm_integral[i] -= Ek_new - Ek;
1272 fprintf(debug,"TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1273 i,Ek_ref,Ek,Ek_new,ekind->tcstat[i].lambda);
1278 ekind->tcstat[i].lambda = 1.0;
1283 real vrescale_energy(t_grpopts *opts,double therm_integral[])
1289 for(i=0; i<opts->ngtc; i++) {
1290 ener += therm_integral[i];
1296 void rescale_velocities(gmx_ekindata_t *ekind,t_mdatoms *mdatoms,
1297 int start,int end,rvec v[])
1300 t_grp_tcstat *tcstat;
1301 unsigned short *cACC,*cTC;
1306 tcstat = ekind->tcstat;
1311 gstat = ekind->grpstat;
1312 cACC = mdatoms->cACC;
1316 for(n=start; n<end; n++)
1326 /* Only scale the velocity component relative to the COM velocity */
1327 rvec_sub(v[n],gstat[ga].u,vrel);
1328 lg = tcstat[gt].lambda;
1329 for(d=0; d<DIM; d++)
1331 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1338 for(n=start; n<end; n++)
1344 lg = tcstat[gt].lambda;
1345 for(d=0; d<DIM; d++)
1354 /* set target temperatures if we are annealing */
1355 void update_annealing_target_temp(t_grpopts *opts,real t)
1358 real pert,thist=0,x;
1360 for(i=0;i<opts->ngtc;i++) {
1361 npoints = opts->anneal_npoints[i];
1362 switch (opts->annealing[i]) {
1366 /* calculate time modulo the period */
1367 pert = opts->anneal_time[i][npoints-1];
1369 thist = t - n*pert; /* modulo time */
1370 /* Make sure rounding didn't get us outside the interval */
1371 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1378 gmx_fatal(FARGS,"Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",i,opts->ngtc,npoints);
1380 /* We are doing annealing for this group if we got here,
1381 * and we have the (relative) time as thist.
1382 * calculate target temp */
1384 while ((j < npoints-1) && (thist>(opts->anneal_time[i][j+1])))
1386 if (j < npoints-1) {
1387 /* Found our position between points j and j+1.
1388 * Interpolate: x is the amount from j+1, (1-x) from point j
1389 * First treat possible jumps in temperature as a special case.
1391 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1392 opts->ref_t[i]=opts->anneal_temp[i][j+1];
1394 x = ((thist-opts->anneal_time[i][j])/
1395 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1396 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1400 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];