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
46 #include "gmx_fatal.h"
49 #include "gmx_random.h"
53 #define NTROTTERPARTS 3
55 /* Suzuki-Yoshida Constants, for n=3 and n=5, for symplectic integration */
57 /* for n=3, w0 = w2 = 1/(2-2^-(1/3)), w1 = 1-2*w0 */
58 /* for n=5, w0 = w1 = w3 = w4 = 1/(4-4^-(1/3)), w1 = 1-4*w0 */
60 #define MAX_SUZUKI_YOSHIDA_NUM 5
61 #define SUZUKI_YOSHIDA_NUM 5
63 static const double sy_const_1[] = { 1. };
64 static const double sy_const_3[] = { 0.828981543588751,-0.657963087177502,0.828981543588751 };
65 static const double sy_const_5[] = { 0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065 };
67 static const double* sy_const[] = {
77 static const double sy_const[MAX_SUZUKI_YOSHIDA_NUM+1][MAX_SUZUKI_YOSHIDA_NUM+1] = {
81 {0.828981543588751,-0.657963087177502,0.828981543588751},
83 {0.2967324292201065,0.2967324292201065,-0.186929716880426,0.2967324292201065,0.2967324292201065}
86 /* these integration routines are only referenced inside this file */
87 static void NHC_trotter(t_grpopts *opts,int nvar, gmx_ekindata_t *ekind,real dtfull,
88 double xi[],double vxi[], double scalefac[], real *veta, t_extmass *MassQ, gmx_bool bEkinAveVel)
91 /* general routine for both barostat and thermostat nose hoover chains */
94 double Ekin,Efac,reft,kT,nd;
101 int mstepsi, mstepsj;
102 int ns = SUZUKI_YOSHIDA_NUM; /* set the degree of integration in the types/state.h file */
103 int nh = opts->nhchainlength;
106 mstepsi = mstepsj = ns;
108 /* if scalefac is NULL, we are doing the NHC of the barostat */
111 if (scalefac == NULL) {
115 for (i=0; i<nvar; i++)
118 /* make it easier to iterate by selecting
119 out the sub-array that corresponds to this T group */
124 iQinv = &(MassQ->QPinv[i*nh]);
125 nd = 1; /* THIS WILL CHANGE IF NOT ISOTROPIC */
126 reft = max(0.0,opts->ref_t[0]);
127 Ekin = sqr(*veta)/MassQ->Winv;
129 iQinv = &(MassQ->Qinv[i*nh]);
130 tcstat = &ekind->tcstat[i];
132 reft = max(0.0,opts->ref_t[i]);
135 Ekin = 2*trace(tcstat->ekinf)*tcstat->ekinscalef_nhc;
137 Ekin = 2*trace(tcstat->ekinh)*tcstat->ekinscaleh_nhc;
142 for(mi=0;mi<mstepsi;mi++)
144 for(mj=0;mj<mstepsj;mj++)
146 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
147 dt = sy_const[ns][mj] * dtfull / mstepsi;
149 /* compute the thermal forces */
150 GQ[0] = iQinv[0]*(Ekin - nd*kT);
154 if (iQinv[j+1] > 0) {
155 /* we actually don't need to update here if we save the
156 state of the GQ, but it's easier to just recompute*/
157 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
163 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
166 Efac = exp(-0.125*dt*ivxi[j]);
167 ivxi[j-1] = Efac*(ivxi[j-1]*Efac + 0.25*dt*GQ[j-1]);
170 Efac = exp(-0.5*dt*ivxi[0]);
178 /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
179 able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
180 think about this a bit more . . . */
182 GQ[0] = iQinv[0]*(Ekin - nd*kT);
184 /* update thermostat positions */
187 ixi[j] += 0.5*dt*ivxi[j];
192 Efac = exp(-0.125*dt*ivxi[j+1]);
193 ivxi[j] = Efac*(ivxi[j]*Efac + 0.25*dt*GQ[j]);
194 if (iQinv[j+1] > 0) {
195 GQ[j+1] = iQinv[j+1]*((sqr(ivxi[j])/iQinv[j])-kT);
200 ivxi[nh-1] += 0.25*dt*GQ[nh-1];
207 static void boxv_trotter(t_inputrec *ir, real *veta, real dt, tensor box,
208 gmx_ekindata_t *ekind, tensor vir, real pcorr, real ecorr, t_extmass *MassQ)
215 tensor Winvm,ekinmod,localpres;
217 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
218 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
221 if (ir->epct==epctSEMIISOTROPIC)
230 /* eta is in pure units. veta is in units of ps^-1. GW is in
231 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
232 taken to use only RATIOS of eta in updating the volume. */
234 /* we take the partial pressure tensors, modify the
235 kinetic energy tensor, and recovert to pressure */
237 if (ir->opts.nrdf[0]==0)
239 gmx_fatal(FARGS,"Barostat is coupled to a T-group with no degrees of freedom\n");
241 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
242 alpha = 1.0 + DIM/((double)ir->opts.nrdf[0]);
243 alpha *= ekind->tcstat[0].ekinscalef_nhc;
244 msmul(ekind->ekin,alpha,ekinmod);
246 pscal = calc_pres(ir->ePBC,nwall,box,ekinmod,vir,localpres) + pcorr;
249 GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
255 * This file implements temperature and pressure coupling algorithms:
256 * For now only the Weak coupling and the modified weak coupling.
258 * Furthermore computation of pressure and temperature is done here
262 real calc_pres(int ePBC,int nwall,matrix box,tensor ekin,tensor vir,
268 if (ePBC==epbcNONE || (ePBC==epbcXY && nwall!=2))
271 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
272 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
276 fac=PRESFAC*2.0/det(box);
277 for(n=0; (n<DIM); n++)
278 for(m=0; (m<DIM); m++)
279 pres[n][m] = (ekin[n][m] - vir[n][m])*fac;
282 pr_rvecs(debug,0,"PC: pres",pres,DIM);
283 pr_rvecs(debug,0,"PC: ekin",ekin,DIM);
284 pr_rvecs(debug,0,"PC: vir ",vir, DIM);
285 pr_rvecs(debug,0,"PC: box ",box, DIM);
288 return trace(pres)/DIM;
291 real calc_temp(real ekin,real nrdf)
294 return (2.0*ekin)/(nrdf*BOLTZ);
299 void parrinellorahman_pcoupl(FILE *fplog,gmx_large_int_t step,
300 t_inputrec *ir,real dt,tensor pres,
301 tensor box,tensor box_rel,tensor boxv,
302 tensor M,matrix mu,gmx_bool bFirstStep)
304 /* This doesn't do any coordinate updating. It just
305 * integrates the box vector equations from the calculated
306 * acceleration due to pressure difference. We also compute
307 * the tensor M which is used in update to couple the particle
308 * coordinates to the box vectors.
310 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
313 * M_nk = (h') * (h' * h + h' h) * h
315 * with the dots denoting time derivatives and h is the transformation from
316 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
317 * This also goes for the pressure and M tensors - they are transposed relative
318 * to ours. Our equation thus becomes:
321 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
323 * where b is the gromacs box matrix.
324 * Our box accelerations are given by
326 * b = vol/W inv(box') * (P-ref_P) (=h')
331 real vol=box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
332 real atot,arel,change,maxchange,xy_pressure;
333 tensor invbox,pdiff,t1,t2;
337 m_inv_ur0(box,invbox);
340 /* Note that PRESFAC does not occur here.
341 * The pressure and compressibility always occur as a product,
342 * therefore the pressure unit drops out.
344 maxl=max(box[XX][XX],box[YY][YY]);
345 maxl=max(maxl,box[ZZ][ZZ]);
349 (4*M_PI*M_PI*ir->compress[d][n])/(3*ir->tau_p*ir->tau_p*maxl);
351 m_sub(pres,ir->ref_p,pdiff);
353 if(ir->epct==epctSURFACETENSION) {
354 /* Unlike Berendsen coupling it might not be trivial to include a z
355 * pressure correction here? On the other hand we don't scale the
356 * box momentarily, but change accelerations, so it might not be crucial.
358 xy_pressure=0.5*(pres[XX][XX]+pres[YY][YY]);
360 pdiff[d][d]=(xy_pressure-(pres[ZZ][ZZ]-ir->ref_p[d][d]/box[d][d]));
363 tmmul(invbox,pdiff,t1);
364 /* Move the off-diagonal elements of the 'force' to one side to ensure
365 * that we obey the box constraints.
369 t1[d][n] += t1[n][d];
375 case epctANISOTROPIC:
378 t1[d][n] *= winv[d][n]*vol;
381 /* calculate total volume acceleration */
382 atot=box[XX][XX]*box[YY][YY]*t1[ZZ][ZZ]+
383 box[XX][XX]*t1[YY][YY]*box[ZZ][ZZ]+
384 t1[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
386 /* set all RELATIVE box accelerations equal, and maintain total V
390 t1[d][n] = winv[0][0]*vol*arel*box[d][n];
392 case epctSEMIISOTROPIC:
393 case epctSURFACETENSION:
394 /* Note the correction to pdiff above for surftens. coupling */
396 /* calculate total XY volume acceleration */
397 atot=box[XX][XX]*t1[YY][YY]+t1[XX][XX]*box[YY][YY];
398 arel=atot/(2*box[XX][XX]*box[YY][YY]);
399 /* set RELATIVE XY box accelerations equal, and maintain total V
400 * change speed. Dont change the third box vector accelerations */
403 t1[d][n] = winv[d][n]*vol*arel*box[d][n];
405 t1[ZZ][n] *= winv[d][n]*vol;
408 gmx_fatal(FARGS,"Parrinello-Rahman pressure coupling type %s "
409 "not supported yet\n",EPCOUPLTYPETYPE(ir->epct));
416 boxv[d][n] += dt*t1[d][n];
418 /* We do NOT update the box vectors themselves here, since
419 * we need them for shifting later. It is instead done last
420 * in the update() routine.
423 /* Calculate the change relative to diagonal elements-
424 since it's perfectly ok for the off-diagonal ones to
425 be zero it doesn't make sense to check the change relative
429 change=fabs(dt*boxv[d][n]/box[d][d]);
431 if (change>maxchange)
435 if (maxchange > 0.01 && fplog) {
438 "\nStep %s Warning: Pressure scaling more than 1%%. "
439 "This may mean your system\n is not yet equilibrated. "
440 "Use of Parrinello-Rahman pressure coupling during\n"
441 "equilibration can lead to simulation instability, "
442 "and is discouraged.\n",
443 gmx_step_str(step,buf));
447 preserve_box_shape(ir,box_rel,boxv);
449 mtmul(boxv,box,t1); /* t1=boxv * b' */
453 /* Determine the scaling matrix mu for the coordinates */
456 t1[d][n] = box[d][n] + dt*boxv[d][n];
457 preserve_box_shape(ir,box_rel,t1);
458 /* t1 is the box at t+dt, determine mu as the relative change */
459 mmul_ur0(invbox,t1,mu);
462 void berendsen_pcoupl(FILE *fplog,gmx_large_int_t step,
463 t_inputrec *ir,real dt, tensor pres,matrix box,
467 real scalar_pressure, xy_pressure, p_corr_z;
468 char *ptr,buf[STRLEN];
471 * Calculate the scaling matrix mu
475 for(d=0; d<DIM; d++) {
476 scalar_pressure += pres[d][d]/DIM;
478 xy_pressure += pres[d][d]/(DIM-1);
480 /* Pressure is now in bar, everywhere. */
481 #define factor(d,m) (ir->compress[d][m]*dt/ir->tau_p)
483 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
484 * necessary for triclinic scaling
491 mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
494 case epctSEMIISOTROPIC:
496 mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
498 1.0 - factor(ZZ,ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
500 case epctANISOTROPIC:
503 mu[d][n] = (d==n ? 1.0 : 0.0)
504 -factor(d,n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
506 case epctSURFACETENSION:
507 /* ir->ref_p[0/1] is the reference surface-tension times *
508 * the number of surfaces */
509 if (ir->compress[ZZ][ZZ])
510 p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
512 /* when the compressibity is zero, set the pressure correction *
513 * in the z-direction to zero to get the correct surface tension */
515 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
516 for(d=0; d<DIM-1; d++)
517 mu[d][d] = 1.0 + factor(d,d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
518 - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
521 gmx_fatal(FARGS,"Berendsen pressure coupling type %s not supported yet\n",
522 EPCOUPLTYPETYPE(ir->epct));
525 /* To fullfill the orientation restrictions on triclinic boxes
526 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
527 * the other elements of mu to first order.
529 mu[YY][XX] += mu[XX][YY];
530 mu[ZZ][XX] += mu[XX][ZZ];
531 mu[ZZ][YY] += mu[YY][ZZ];
537 pr_rvecs(debug,0,"PC: pres ",pres,3);
538 pr_rvecs(debug,0,"PC: mu ",mu,3);
541 if (mu[XX][XX]<0.99 || mu[XX][XX]>1.01 ||
542 mu[YY][YY]<0.99 || mu[YY][YY]>1.01 ||
543 mu[ZZ][ZZ]<0.99 || mu[ZZ][ZZ]>1.01) {
545 sprintf(buf,"\nStep %s Warning: pressure scaling more than 1%%, "
547 gmx_step_str(step,buf2),mu[XX][XX],mu[YY][YY],mu[ZZ][ZZ]);
549 fprintf(fplog,"%s",buf);
550 fprintf(stderr,"%s",buf);
554 void berendsen_pscale(t_inputrec *ir,matrix mu,
555 matrix box,matrix box_rel,
556 int start,int nr_atoms,
557 rvec x[],unsigned short cFREEZE[],
560 ivec *nFreeze=ir->opts.nFreeze;
563 /* Scale the positions */
564 for (n=start; n<start+nr_atoms; n++) {
569 x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
571 x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
573 x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
575 /* compute final boxlengths */
576 for (d=0; d<DIM; d++) {
577 box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
578 box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
579 box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
582 preserve_box_shape(ir,box_rel,box);
584 /* (un)shifting should NOT be done after this,
585 * since the box vectors might have changed
587 inc_nrnb(nrnb,eNR_PCOUPL,nr_atoms);
590 void berendsen_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt)
598 for(i=0; (i<opts->ngtc); i++)
602 T = ekind->tcstat[i].T;
606 T = ekind->tcstat[i].Th;
609 if ((opts->tau_t[i] > 0) && (T > 0.0)) {
611 reft = max(0.0,opts->ref_t[i]);
612 lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
613 ekind->tcstat[i].lambda = max(min(lll,1.25),0.8);
616 ekind->tcstat[i].lambda = 1.0;
620 fprintf(debug,"TC: group %d: T: %g, Lambda: %g\n",
621 i,T,ekind->tcstat[i].lambda);
625 void nosehoover_tcoupl(t_grpopts *opts,gmx_ekindata_t *ekind,real dt,
626 double xi[],double vxi[], t_extmass *MassQ)
631 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
633 for(i=0; (i<opts->ngtc); i++) {
634 reft = max(0.0,opts->ref_t[i]);
636 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
637 xi[i] += dt*(oldvxi + vxi[i])*0.5;
641 t_state *init_bufstate(const t_state *template_state)
644 int nc = template_state->nhchainlength;
646 snew(state->nosehoover_xi,nc*template_state->ngtc);
647 snew(state->nosehoover_vxi,nc*template_state->ngtc);
648 snew(state->therm_integral,template_state->ngtc);
649 snew(state->nhpres_xi,nc*template_state->nnhpres);
650 snew(state->nhpres_vxi,nc*template_state->nnhpres);
655 void destroy_bufstate(t_state *state)
659 sfree(state->nosehoover_xi);
660 sfree(state->nosehoover_vxi);
661 sfree(state->therm_integral);
662 sfree(state->nhpres_xi);
663 sfree(state->nhpres_vxi);
667 void trotter_update(t_inputrec *ir,gmx_large_int_t step, gmx_ekindata_t *ekind,
668 gmx_enerdata_t *enerd, t_state *state,
669 tensor vir, t_mdatoms *md,
670 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
673 int n,i,j,d,ntgrp,ngtc,gc=0;
674 t_grp_tcstat *tcstat;
676 gmx_large_int_t step_eff;
677 real ecorr,pcorr,dvdlcorr;
678 real bmass,qmass,reft,kT,dt,nd;
679 tensor dumpres,dumvir;
680 double *scalefac,dtc;
682 rvec sumv={0,0,0},consk;
685 if (trotter_seqno <= ettTSEQ2)
687 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
688 is actually the last half step from the previous step. Thus the first half step
689 actually corresponds to the n-1 step*/
695 bCouple = (ir->nsttcouple == 1 ||
696 do_per_step(step_eff+ir->nsttcouple,ir->nsttcouple));
698 trotter_seq = trotter_seqlist[trotter_seqno];
700 /* signal we are returning if nothing is going to be done in this routine */
701 if ((trotter_seq[0] == etrtSKIPALL) || !(bCouple))
706 dtc = ir->nsttcouple*ir->delta_t;
707 opts = &(ir->opts); /* just for ease of referencing */
709 snew(scalefac,opts->ngtc);
714 /* execute the series of trotter updates specified in the trotterpart array */
716 for (i=0;i<NTROTTERPARTS;i++){
717 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
718 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
727 switch (trotter_seq[i])
731 boxv_trotter(ir,&(state->veta),dt,state->box,ekind,vir,
732 enerd->term[F_PDISPCORR],enerd->term[F_DISPCORR],MassQ);
736 NHC_trotter(opts,state->nnhpres,ekind,dt,state->nhpres_xi,
737 state->nhpres_vxi,NULL,&(state->veta),MassQ,FALSE);
741 NHC_trotter(opts,opts->ngtc,ekind,dt,state->nosehoover_xi,
742 state->nosehoover_vxi,scalefac,NULL,MassQ,(ir->eI==eiVV));
743 /* need to rescale the kinetic energies and velocities here. Could
744 scale the velocities later, but we need them scaled in order to
745 produce the correct outputs, so we'll scale them here. */
747 for (i=0; i<ngtc;i++)
749 tcstat = &ekind->tcstat[i];
750 tcstat->vscale_nhc = scalefac[i];
751 tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
752 tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
754 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
755 /* but do we actually need the total? */
757 /* modify the velocities as well */
758 for (n=md->start;n<md->start+md->homenr;n++)
766 state->v[n][d] *= scalefac[gc];
773 sumv[d] += (state->v[n][d])/md->invmass[n];
782 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
790 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
792 fprintf(debug,"Conserved kappa: %15.8f %15.8f %15.8f\n",consk[0],consk[1],consk[2]);
799 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
801 int n,i,j,d,ntgrp,ngtc,nnhpres,nh,gc=0;
802 t_grp_tcstat *tcstat;
804 real ecorr,pcorr,dvdlcorr;
805 real bmass,qmass,reft,kT,dt,ndj,nd;
806 tensor dumpres,dumvir;
809 opts = &(ir->opts); /* just for ease of referencing */
811 nnhpres = state->nnhpres;
812 nh = state->nhchainlength;
814 if (ir->eI == eiMD) {
815 snew(MassQ->Qinv,ngtc);
816 for(i=0; (i<ngtc); i++)
818 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
820 MassQ->Qinv[i]=1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
828 else if (EI_VV(ir->eI))
830 /* Set pressure variables */
832 if (state->vol0 == 0)
834 state->vol0 = det(state->box); /* because we start by defining a fixed compressibility,
835 we need the volume at this compressibility to solve the problem */
838 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
839 /* Investigate this more -- is this the right mass to make it? */
840 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
841 /* An alternate mass definition, from Tuckerman et al. */
842 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
847 MassQ->Winvm[d][n]= PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
848 /* not clear this is correct yet for the anisotropic case*/
851 /* Allocate space for thermostat variables */
852 snew(MassQ->Qinv,ngtc*nh);
854 /* now, set temperature variables */
855 for(i=0; i<ngtc; i++)
857 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
859 reft = max(0.0,opts->ref_t[i]);
872 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
880 MassQ->Qinv[i*nh+j] = 0.0;
886 /* first, initialize clear all the trotter calls */
887 snew(trotter_seq,ettTSEQMAX);
888 for (i=0;i<ettTSEQMAX;i++)
890 snew(trotter_seq[i],NTROTTERPARTS);
891 for (j=0;j<NTROTTERPARTS;j++) {
892 trotter_seq[i][j] = etrtNONE;
894 trotter_seq[i][0] = etrtSKIPALL;
899 /* no trotter calls, so we never use the values in the array.
900 * We access them (so we need to define them, but ignore
906 /* compute the kinetic energy by using the half step velocities or
907 * the kinetic energies, depending on the order of the trotter calls */
911 if (IR_NPT_TROTTER(ir))
913 /* This is the complicated version - there are 4 possible calls, depending on ordering.
914 We start with the initial one. */
915 /* first, a round that estimates veta. */
916 trotter_seq[0][0] = etrtBAROV;
918 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
920 /* The first half trotter update */
921 trotter_seq[2][0] = etrtBAROV;
922 trotter_seq[2][1] = etrtNHC;
923 trotter_seq[2][2] = etrtBARONHC;
925 /* The second half trotter update */
926 trotter_seq[3][0] = etrtBARONHC;
927 trotter_seq[3][1] = etrtNHC;
928 trotter_seq[3][2] = etrtBAROV;
930 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
935 if (IR_NVT_TROTTER(ir))
937 /* This is the easy version - there are only two calls, both the same.
938 Otherwise, even easier -- no calls */
939 trotter_seq[2][0] = etrtNHC;
940 trotter_seq[3][0] = etrtNHC;
943 } else if (ir->eI==eiVVAK) {
944 if (IR_NPT_TROTTER(ir))
946 /* This is the complicated version - there are 4 possible calls, depending on ordering.
947 We start with the initial one. */
948 /* first, a round that estimates veta. */
949 trotter_seq[0][0] = etrtBAROV;
951 /* The first half trotter update, part 1 -- double update, because it commutes */
952 trotter_seq[1][0] = etrtNHC;
954 /* The first half trotter update, part 2 */
955 trotter_seq[2][0] = etrtBAROV;
956 trotter_seq[2][1] = etrtBARONHC;
958 /* The second half trotter update, part 1 */
959 trotter_seq[3][0] = etrtBARONHC;
960 trotter_seq[3][1] = etrtBAROV;
962 /* The second half trotter update -- blank for now */
963 trotter_seq[4][0] = etrtNHC;
967 if (IR_NVT_TROTTER(ir))
969 /* This is the easy version - there is only one call, both the same.
970 Otherwise, even easier -- no calls */
971 trotter_seq[1][0] = etrtNHC;
972 trotter_seq[4][0] = etrtNHC;
981 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
984 snew(MassQ->QPinv,nnhpres*opts->nhchainlength);
986 /* barostat temperature */
987 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
989 reft = max(0.0,opts->ref_t[0]);
991 for (i=0;i<nnhpres;i++) {
1001 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1007 for (i=0;i<nnhpres;i++) {
1010 MassQ->QPinv[i*nh+j]=0.0;
1017 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1019 int i,j,nd,ndj,bmass,qmass,ngtcall;
1020 real ener_npt,reft,eta,kT,tau;
1024 int nh = state->nhchainlength;
1028 /* now we compute the contribution of the pressure to the conserved quantity*/
1030 if (ir->epc==epcMTTK)
1032 /* find the volume, and the kinetic energy of the volume */
1037 /* contribution from the pressure momenenta */
1038 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1040 /* contribution from the PV term */
1041 vol = det(state->box);
1042 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1045 case epctANISOTROPIC:
1049 case epctSURFACETENSION:
1052 case epctSEMIISOTROPIC:
1060 if (IR_NPT_TROTTER(ir))
1062 /* add the energy from the barostat thermostat chain */
1063 for (i=0;i<state->nnhpres;i++) {
1065 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1066 ivxi = &state->nhpres_vxi[i*nh];
1067 ixi = &state->nhpres_xi[i*nh];
1068 iQinv = &(MassQ->QPinv[i*nh]);
1069 reft = max(ir->opts.ref_t[0],0); /* using 'System' temperature */
1076 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1077 /* contribution from the thermal variable of the NH chain */
1078 ener_npt += ixi[j]*kT;
1082 fprintf(debug,"P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f",i,j,ivxi[j],ixi[j]);
1090 for(i=0; i<ir->opts.ngtc; i++)
1092 ixi = &state->nosehoover_xi[i*nh];
1093 ivxi = &state->nosehoover_vxi[i*nh];
1094 iQinv = &(MassQ->Qinv[i*nh]);
1096 nd = ir->opts.nrdf[i];
1097 reft = max(ir->opts.ref_t[i],0);
1102 if (IR_NVT_TROTTER(ir))
1104 /* contribution from the thermal momenta of the NH chain */
1108 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1109 /* contribution from the thermal variable of the NH chain */
1117 ener_npt += ndj*ixi[j]*kT;
1121 else /* Other non Trotter temperature NH control -- no chains yet. */
1123 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1124 ener_npt += nd*ixi[0]*kT;
1132 static real vrescale_gamdev(int ia, gmx_rng_t rng)
1133 /* Gamma distribution, adapted from numerical recipes */
1136 real am,e,s,v1,v2,x,y;
1143 for(j=1; j<=ia; j++)
1145 x *= gmx_rng_uniform_real(rng);
1159 v1 = gmx_rng_uniform_real(rng);
1160 v2 = 2.0*gmx_rng_uniform_real(rng)-1.0;
1162 while (v1*v1 + v2*v2 > 1.0 ||
1163 v1*v1*GMX_REAL_MAX < 3.0*ia);
1164 /* The last check above ensures that both x (3.0 > 2.0 in s)
1165 * and the pre-factor for e do not go out of range.
1169 s = sqrt(2.0*am + 1.0);
1173 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1175 while (gmx_rng_uniform_real(rng) > e);
1181 static real vrescale_sumnoises(int nn,gmx_rng_t rng)
1184 * Returns the sum of n independent gaussian noises squared
1185 * (i.e. equivalent to summing the square of the return values
1186 * of nn calls to gmx_rng_gaussian_real).xs
1192 } else if (nn == 1) {
1193 rr = gmx_rng_gaussian_real(rng);
1195 } else if (nn % 2 == 0) {
1196 return 2.0*vrescale_gamdev(nn/2,rng);
1198 rr = gmx_rng_gaussian_real(rng);
1199 return 2.0*vrescale_gamdev((nn-1)/2,rng) + rr*rr;
1203 static real vrescale_resamplekin(real kk,real sigma, int ndeg, real taut,
1207 * Generates a new value for the kinetic energy,
1208 * according to Bussi et al JCP (2007), Eq. (A7)
1209 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1210 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1211 * ndeg: number of degrees of freedom of the atoms to be thermalized
1212 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1217 factor = exp(-1.0/taut);
1221 rr = gmx_rng_gaussian_real(rng);
1224 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1,rng) + rr*rr)/ndeg - kk) +
1225 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1228 void vrescale_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt,
1229 double therm_integral[],gmx_rng_t rng)
1233 real Ek,Ek_ref1,Ek_ref,Ek_new;
1237 for(i=0; (i<opts->ngtc); i++)
1241 Ek = trace(ekind->tcstat[i].ekinf);
1245 Ek = trace(ekind->tcstat[i].ekinh);
1248 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1250 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1251 Ek_ref = Ek_ref1*opts->nrdf[i];
1253 Ek_new = vrescale_resamplekin(Ek,Ek_ref,opts->nrdf[i],
1254 opts->tau_t[i]/dt,rng);
1256 /* Analytically Ek_new>=0, but we check for rounding errors */
1259 ekind->tcstat[i].lambda = 0.0;
1263 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1266 therm_integral[i] -= Ek_new - Ek;
1270 fprintf(debug,"TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1271 i,Ek_ref,Ek,Ek_new,ekind->tcstat[i].lambda);
1276 ekind->tcstat[i].lambda = 1.0;
1281 real vrescale_energy(t_grpopts *opts,double therm_integral[])
1287 for(i=0; i<opts->ngtc; i++) {
1288 ener += therm_integral[i];
1294 void rescale_velocities(gmx_ekindata_t *ekind,t_mdatoms *mdatoms,
1295 int start,int end,rvec v[])
1298 t_grp_tcstat *tcstat;
1299 unsigned short *cACC,*cTC;
1304 tcstat = ekind->tcstat;
1309 gstat = ekind->grpstat;
1310 cACC = mdatoms->cACC;
1314 for(n=start; n<end; n++)
1324 /* Only scale the velocity component relative to the COM velocity */
1325 rvec_sub(v[n],gstat[ga].u,vrel);
1326 lg = tcstat[gt].lambda;
1327 for(d=0; d<DIM; d++)
1329 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1336 for(n=start; n<end; n++)
1342 lg = tcstat[gt].lambda;
1343 for(d=0; d<DIM; d++)
1352 /* set target temperatures if we are annealing */
1353 void update_annealing_target_temp(t_grpopts *opts,real t)
1356 real pert,thist=0,x;
1358 for(i=0;i<opts->ngtc;i++) {
1359 npoints = opts->anneal_npoints[i];
1360 switch (opts->annealing[i]) {
1364 /* calculate time modulo the period */
1365 pert = opts->anneal_time[i][npoints-1];
1367 thist = t - n*pert; /* modulo time */
1368 /* Make sure rounding didn't get us outside the interval */
1369 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1376 gmx_fatal(FARGS,"Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",i,opts->ngtc,npoints);
1378 /* We are doing annealing for this group if we got here,
1379 * and we have the (relative) time as thist.
1380 * calculate target temp */
1382 while ((j < npoints-1) && (thist>(opts->anneal_time[i][j+1])))
1384 if (j < npoints-1) {
1385 /* Found our position between points j and j+1.
1386 * Interpolate: x is the amount from j+1, (1-x) from point j
1387 * First treat possible jumps in temperature as a special case.
1389 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1390 opts->ref_t[i]=opts->anneal_temp[i][j+1];
1392 x = ((thist-opts->anneal_time[i][j])/
1393 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1394 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1398 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];