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 /* for now, we use Elr = 0, because if you want to get it right, you
247 really should be using PME. Maybe print a warning? */
249 pscal = calc_pres(ir->ePBC,nwall,box,ekinmod,vir,localpres,0.0) + pcorr;
252 GW = (vol*(MassQ->Winv/PRESFAC))*(DIM*pscal - trace(ir->ref_p)); /* W is in ps^2 * bar * nm^3 */
258 * This file implements temperature and pressure coupling algorithms:
259 * For now only the Weak coupling and the modified weak coupling.
261 * Furthermore computation of pressure and temperature is done here
265 real calc_pres(int ePBC,int nwall,matrix box,tensor ekin,tensor vir,
266 tensor pres,real Elr)
271 if (ePBC==epbcNONE || (ePBC==epbcXY && nwall!=2))
274 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
275 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
279 /* Long range correction for periodic systems, see
281 * divide by 6 because it is multiplied by fac later on.
282 * If Elr = 0, no correction is made.
285 /* This formula should not be used with Ewald or PME,
286 * where the full long-range virial is calculated. EL 990823
290 fac=PRESFAC*2.0/det(box);
291 for(n=0; (n<DIM); n++)
292 for(m=0; (m<DIM); m++)
293 pres[n][m]=(ekin[n][m]-vir[n][m]+Plr)*fac;
296 pr_rvecs(debug,0,"PC: pres",pres,DIM);
297 pr_rvecs(debug,0,"PC: ekin",ekin,DIM);
298 pr_rvecs(debug,0,"PC: vir ",vir, DIM);
299 pr_rvecs(debug,0,"PC: box ",box, DIM);
302 return trace(pres)/DIM;
305 real calc_temp(real ekin,real nrdf)
308 return (2.0*ekin)/(nrdf*BOLTZ);
313 void parrinellorahman_pcoupl(FILE *fplog,gmx_large_int_t step,
314 t_inputrec *ir,real dt,tensor pres,
315 tensor box,tensor box_rel,tensor boxv,
316 tensor M,matrix mu,gmx_bool bFirstStep)
318 /* This doesn't do any coordinate updating. It just
319 * integrates the box vector equations from the calculated
320 * acceleration due to pressure difference. We also compute
321 * the tensor M which is used in update to couple the particle
322 * coordinates to the box vectors.
324 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
327 * M_nk = (h') * (h' * h + h' h) * h
329 * with the dots denoting time derivatives and h is the transformation from
330 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
331 * This also goes for the pressure and M tensors - they are transposed relative
332 * to ours. Our equation thus becomes:
335 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
337 * where b is the gromacs box matrix.
338 * Our box accelerations are given by
340 * b = vol/W inv(box') * (P-ref_P) (=h')
345 real vol=box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
346 real atot,arel,change,maxchange,xy_pressure;
347 tensor invbox,pdiff,t1,t2;
351 m_inv_ur0(box,invbox);
354 /* Note that PRESFAC does not occur here.
355 * The pressure and compressibility always occur as a product,
356 * therefore the pressure unit drops out.
358 maxl=max(box[XX][XX],box[YY][YY]);
359 maxl=max(maxl,box[ZZ][ZZ]);
363 (4*M_PI*M_PI*ir->compress[d][n])/(3*ir->tau_p*ir->tau_p*maxl);
365 m_sub(pres,ir->ref_p,pdiff);
367 if(ir->epct==epctSURFACETENSION) {
368 /* Unlike Berendsen coupling it might not be trivial to include a z
369 * pressure correction here? On the other hand we don't scale the
370 * box momentarily, but change accelerations, so it might not be crucial.
372 xy_pressure=0.5*(pres[XX][XX]+pres[YY][YY]);
374 pdiff[d][d]=(xy_pressure-(pres[ZZ][ZZ]-ir->ref_p[d][d]/box[d][d]));
377 tmmul(invbox,pdiff,t1);
378 /* Move the off-diagonal elements of the 'force' to one side to ensure
379 * that we obey the box constraints.
383 t1[d][n] += t1[n][d];
389 case epctANISOTROPIC:
392 t1[d][n] *= winv[d][n]*vol;
395 /* calculate total volume acceleration */
396 atot=box[XX][XX]*box[YY][YY]*t1[ZZ][ZZ]+
397 box[XX][XX]*t1[YY][YY]*box[ZZ][ZZ]+
398 t1[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
400 /* set all RELATIVE box accelerations equal, and maintain total V
404 t1[d][n] = winv[0][0]*vol*arel*box[d][n];
406 case epctSEMIISOTROPIC:
407 case epctSURFACETENSION:
408 /* Note the correction to pdiff above for surftens. coupling */
410 /* calculate total XY volume acceleration */
411 atot=box[XX][XX]*t1[YY][YY]+t1[XX][XX]*box[YY][YY];
412 arel=atot/(2*box[XX][XX]*box[YY][YY]);
413 /* set RELATIVE XY box accelerations equal, and maintain total V
414 * change speed. Dont change the third box vector accelerations */
417 t1[d][n] = winv[d][n]*vol*arel*box[d][n];
419 t1[ZZ][n] *= winv[d][n]*vol;
422 gmx_fatal(FARGS,"Parrinello-Rahman pressure coupling type %s "
423 "not supported yet\n",EPCOUPLTYPETYPE(ir->epct));
430 boxv[d][n] += dt*t1[d][n];
432 /* We do NOT update the box vectors themselves here, since
433 * we need them for shifting later. It is instead done last
434 * in the update() routine.
437 /* Calculate the change relative to diagonal elements-
438 since it's perfectly ok for the off-diagonal ones to
439 be zero it doesn't make sense to check the change relative
443 change=fabs(dt*boxv[d][n]/box[d][d]);
445 if (change>maxchange)
449 if (maxchange > 0.01 && fplog) {
452 "\nStep %s Warning: Pressure scaling more than 1%%. "
453 "This may mean your system\n is not yet equilibrated. "
454 "Use of Parrinello-Rahman pressure coupling during\n"
455 "equilibration can lead to simulation instability, "
456 "and is discouraged.\n",
457 gmx_step_str(step,buf));
461 preserve_box_shape(ir,box_rel,boxv);
463 mtmul(boxv,box,t1); /* t1=boxv * b' */
467 /* Determine the scaling matrix mu for the coordinates */
470 t1[d][n] = box[d][n] + dt*boxv[d][n];
471 preserve_box_shape(ir,box_rel,t1);
472 /* t1 is the box at t+dt, determine mu as the relative change */
473 mmul_ur0(invbox,t1,mu);
476 void berendsen_pcoupl(FILE *fplog,gmx_large_int_t step,
477 t_inputrec *ir,real dt, tensor pres,matrix box,
481 real scalar_pressure, xy_pressure, p_corr_z;
482 char *ptr,buf[STRLEN];
485 * Calculate the scaling matrix mu
489 for(d=0; d<DIM; d++) {
490 scalar_pressure += pres[d][d]/DIM;
492 xy_pressure += pres[d][d]/(DIM-1);
494 /* Pressure is now in bar, everywhere. */
495 #define factor(d,m) (ir->compress[d][m]*dt/ir->tau_p)
497 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
498 * necessary for triclinic scaling
505 mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
508 case epctSEMIISOTROPIC:
510 mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
512 1.0 - factor(ZZ,ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
514 case epctANISOTROPIC:
517 mu[d][n] = (d==n ? 1.0 : 0.0)
518 -factor(d,n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
520 case epctSURFACETENSION:
521 /* ir->ref_p[0/1] is the reference surface-tension times *
522 * the number of surfaces */
523 if (ir->compress[ZZ][ZZ])
524 p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
526 /* when the compressibity is zero, set the pressure correction *
527 * in the z-direction to zero to get the correct surface tension */
529 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
530 for(d=0; d<DIM-1; d++)
531 mu[d][d] = 1.0 + factor(d,d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
532 - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
535 gmx_fatal(FARGS,"Berendsen pressure coupling type %s not supported yet\n",
536 EPCOUPLTYPETYPE(ir->epct));
539 /* To fullfill the orientation restrictions on triclinic boxes
540 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
541 * the other elements of mu to first order.
543 mu[YY][XX] += mu[XX][YY];
544 mu[ZZ][XX] += mu[XX][ZZ];
545 mu[ZZ][YY] += mu[YY][ZZ];
551 pr_rvecs(debug,0,"PC: pres ",pres,3);
552 pr_rvecs(debug,0,"PC: mu ",mu,3);
555 if (mu[XX][XX]<0.99 || mu[XX][XX]>1.01 ||
556 mu[YY][YY]<0.99 || mu[YY][YY]>1.01 ||
557 mu[ZZ][ZZ]<0.99 || mu[ZZ][ZZ]>1.01) {
559 sprintf(buf,"\nStep %s Warning: pressure scaling more than 1%%, "
561 gmx_step_str(step,buf2),mu[XX][XX],mu[YY][YY],mu[ZZ][ZZ]);
563 fprintf(fplog,"%s",buf);
564 fprintf(stderr,"%s",buf);
568 void berendsen_pscale(t_inputrec *ir,matrix mu,
569 matrix box,matrix box_rel,
570 int start,int nr_atoms,
571 rvec x[],unsigned short cFREEZE[],
574 ivec *nFreeze=ir->opts.nFreeze;
577 /* Scale the positions */
578 for (n=start; n<start+nr_atoms; n++) {
583 x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
585 x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
587 x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
589 /* compute final boxlengths */
590 for (d=0; d<DIM; d++) {
591 box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
592 box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
593 box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
596 preserve_box_shape(ir,box_rel,box);
598 /* (un)shifting should NOT be done after this,
599 * since the box vectors might have changed
601 inc_nrnb(nrnb,eNR_PCOUPL,nr_atoms);
604 void berendsen_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt)
612 for(i=0; (i<opts->ngtc); i++)
616 T = ekind->tcstat[i].T;
620 T = ekind->tcstat[i].Th;
623 if ((opts->tau_t[i] > 0) && (T > 0.0)) {
625 reft = max(0.0,opts->ref_t[i]);
626 lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
627 ekind->tcstat[i].lambda = max(min(lll,1.25),0.8);
630 ekind->tcstat[i].lambda = 1.0;
634 fprintf(debug,"TC: group %d: T: %g, Lambda: %g\n",
635 i,T,ekind->tcstat[i].lambda);
639 void nosehoover_tcoupl(t_grpopts *opts,gmx_ekindata_t *ekind,real dt,
640 double xi[],double vxi[], t_extmass *MassQ)
645 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
647 for(i=0; (i<opts->ngtc); i++) {
648 reft = max(0.0,opts->ref_t[i]);
650 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
651 xi[i] += dt*(oldvxi + vxi[i])*0.5;
655 t_state *init_bufstate(const t_state *template_state)
658 int nc = template_state->nhchainlength;
660 snew(state->nosehoover_xi,nc*template_state->ngtc);
661 snew(state->nosehoover_vxi,nc*template_state->ngtc);
662 snew(state->therm_integral,template_state->ngtc);
663 snew(state->nhpres_xi,nc*template_state->nnhpres);
664 snew(state->nhpres_vxi,nc*template_state->nnhpres);
669 void destroy_bufstate(t_state *state)
673 sfree(state->nosehoover_xi);
674 sfree(state->nosehoover_vxi);
675 sfree(state->therm_integral);
676 sfree(state->nhpres_xi);
677 sfree(state->nhpres_vxi);
681 void trotter_update(t_inputrec *ir,gmx_large_int_t step, gmx_ekindata_t *ekind,
682 gmx_enerdata_t *enerd, t_state *state,
683 tensor vir, t_mdatoms *md,
684 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
687 int n,i,j,d,ntgrp,ngtc,gc=0;
688 t_grp_tcstat *tcstat;
690 gmx_large_int_t step_eff;
691 real ecorr,pcorr,dvdlcorr;
692 real bmass,qmass,reft,kT,dt,nd;
693 tensor dumpres,dumvir;
694 double *scalefac,dtc;
699 if (trotter_seqno <= ettTSEQ2)
701 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
702 is actually the last half step from the previous step. Thus the first half step
703 actually corresponds to the n-1 step*/
709 bCouple = (ir->nsttcouple == 1 ||
710 do_per_step(step_eff+ir->nsttcouple,ir->nsttcouple));
712 trotter_seq = trotter_seqlist[trotter_seqno];
714 /* signal we are returning if nothing is going to be done in this routine */
715 if ((trotter_seq[0] == etrtSKIPALL) || !(bCouple))
720 dtc = ir->nsttcouple*ir->delta_t;
721 opts = &(ir->opts); /* just for ease of referencing */
723 snew(scalefac,opts->ngtc);
728 /* execute the series of trotter updates specified in the trotterpart array */
730 for (i=0;i<NTROTTERPARTS;i++){
731 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
732 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
741 switch (trotter_seq[i])
745 boxv_trotter(ir,&(state->veta),dt,state->box,ekind,vir,
746 enerd->term[F_PDISPCORR],enerd->term[F_DISPCORR],MassQ);
750 NHC_trotter(opts,state->nnhpres,ekind,dt,state->nhpres_xi,
751 state->nhpres_vxi,NULL,&(state->veta),MassQ,FALSE);
755 NHC_trotter(opts,opts->ngtc,ekind,dt,state->nosehoover_xi,
756 state->nosehoover_vxi,scalefac,NULL,MassQ,(ir->eI==eiVV));
757 /* need to rescale the kinetic energies and velocities here. Could
758 scale the velocities later, but we need them scaled in order to
759 produce the correct outputs, so we'll scale them here. */
761 for (i=0; i<ngtc;i++)
763 tcstat = &ekind->tcstat[i];
764 tcstat->vscale_nhc = scalefac[i];
765 tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
766 tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
768 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
769 /* but do we actually need the total? */
771 /* modify the velocities as well */
772 for (n=md->start;n<md->start+md->homenr;n++)
780 state->v[n][d] *= scalefac[gc];
787 sumv[d] += (state->v[n][d])/md->invmass[n];
796 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
804 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
806 fprintf(debug,"Conserved kappa: %15.8f %15.8f %15.8f\n",consk[0],consk[1],consk[2]);
813 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
815 int n,i,j,d,ntgrp,ngtc,nnhpres,nh,gc=0;
816 t_grp_tcstat *tcstat;
818 real ecorr,pcorr,dvdlcorr;
819 real bmass,qmass,reft,kT,dt,ndj,nd;
820 tensor dumpres,dumvir;
823 opts = &(ir->opts); /* just for ease of referencing */
825 nnhpres = state->nnhpres;
826 nh = state->nhchainlength;
828 if (ir->eI == eiMD) {
829 snew(MassQ->Qinv,ngtc);
830 for(i=0; (i<ngtc); i++)
832 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
834 MassQ->Qinv[i]=1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
842 else if (EI_VV(ir->eI))
844 /* Set pressure variables */
846 if (state->vol0 == 0)
848 state->vol0 = det(state->box); /* because we start by defining a fixed compressibility,
849 we need the volume at this compressibility to solve the problem */
852 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
853 /* Investigate this more -- is this the right mass to make it? */
854 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
855 /* An alternate mass definition, from Tuckerman et al. */
856 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
861 MassQ->Winvm[d][n]= PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
862 /* not clear this is correct yet for the anisotropic case*/
865 /* Allocate space for thermostat variables */
866 snew(MassQ->Qinv,ngtc*nh);
868 /* now, set temperature variables */
869 for(i=0; i<ngtc; i++)
871 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
873 reft = max(0.0,opts->ref_t[i]);
886 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
894 MassQ->Qinv[i*nh+j] = 0.0;
900 /* first, initialize clear all the trotter calls */
901 snew(trotter_seq,ettTSEQMAX);
902 for (i=0;i<ettTSEQMAX;i++)
904 snew(trotter_seq[i],NTROTTERPARTS);
905 for (j=0;j<NTROTTERPARTS;j++) {
906 trotter_seq[i][j] = etrtNONE;
908 trotter_seq[i][0] = etrtSKIPALL;
913 /* no trotter calls, so we never use the values in the array.
914 * We access them (so we need to define them, but ignore
920 /* compute the kinetic energy by using the half step velocities or
921 * the kinetic energies, depending on the order of the trotter calls */
925 if (IR_NPT_TROTTER(ir))
927 /* This is the complicated version - there are 4 possible calls, depending on ordering.
928 We start with the initial one. */
929 /* first, a round that estimates veta. */
930 trotter_seq[0][0] = etrtBAROV;
932 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
934 /* The first half trotter update */
935 trotter_seq[2][0] = etrtBAROV;
936 trotter_seq[2][1] = etrtNHC;
937 trotter_seq[2][2] = etrtBARONHC;
939 /* The second half trotter update */
940 trotter_seq[3][0] = etrtBARONHC;
941 trotter_seq[3][1] = etrtNHC;
942 trotter_seq[3][2] = etrtBAROV;
944 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
949 if (IR_NVT_TROTTER(ir))
951 /* This is the easy version - there are only two calls, both the same.
952 Otherwise, even easier -- no calls */
953 trotter_seq[2][0] = etrtNHC;
954 trotter_seq[3][0] = etrtNHC;
957 } else if (ir->eI==eiVVAK) {
958 if (IR_NPT_TROTTER(ir))
960 /* This is the complicated version - there are 4 possible calls, depending on ordering.
961 We start with the initial one. */
962 /* first, a round that estimates veta. */
963 trotter_seq[0][0] = etrtBAROV;
965 /* The first half trotter update, part 1 -- double update, because it commutes */
966 trotter_seq[1][0] = etrtNHC;
968 /* The first half trotter update, part 2 */
969 trotter_seq[2][0] = etrtBAROV;
970 trotter_seq[2][1] = etrtBARONHC;
972 /* The second half trotter update, part 1 */
973 trotter_seq[3][0] = etrtBARONHC;
974 trotter_seq[3][1] = etrtBAROV;
976 /* The second half trotter update -- blank for now */
977 trotter_seq[4][0] = etrtNHC;
981 if (IR_NVT_TROTTER(ir))
983 /* This is the easy version - there is only one call, both the same.
984 Otherwise, even easier -- no calls */
985 trotter_seq[1][0] = etrtNHC;
986 trotter_seq[4][0] = etrtNHC;
995 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
998 snew(MassQ->QPinv,nnhpres*opts->nhchainlength);
1000 /* barostat temperature */
1001 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
1003 reft = max(0.0,opts->ref_t[0]);
1005 for (i=0;i<nnhpres;i++) {
1015 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1021 for (i=0;i<nnhpres;i++) {
1024 MassQ->QPinv[i*nh+j]=0.0;
1031 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1033 int i,j,nd,ndj,bmass,qmass,ngtcall;
1034 real ener_npt,reft,eta,kT,tau;
1038 int nh = state->nhchainlength;
1042 /* now we compute the contribution of the pressure to the conserved quantity*/
1044 if (ir->epc==epcMTTK)
1046 /* find the volume, and the kinetic energy of the volume */
1051 /* contribution from the pressure momenenta */
1052 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1054 /* contribution from the PV term */
1055 vol = det(state->box);
1056 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1059 case epctANISOTROPIC:
1063 case epctSURFACETENSION:
1066 case epctSEMIISOTROPIC:
1074 if (IR_NPT_TROTTER(ir))
1076 /* add the energy from the barostat thermostat chain */
1077 for (i=0;i<state->nnhpres;i++) {
1079 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1080 ivxi = &state->nhpres_vxi[i*nh];
1081 ixi = &state->nhpres_xi[i*nh];
1082 iQinv = &(MassQ->QPinv[i*nh]);
1083 reft = max(ir->opts.ref_t[0],0); /* using 'System' temperature */
1090 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1091 /* contribution from the thermal variable of the NH chain */
1092 ener_npt += ixi[j]*kT;
1096 fprintf(debug,"P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f",i,j,ivxi[j],ixi[j]);
1104 for(i=0; i<ir->opts.ngtc; i++)
1106 ixi = &state->nosehoover_xi[i*nh];
1107 ivxi = &state->nosehoover_vxi[i*nh];
1108 iQinv = &(MassQ->Qinv[i*nh]);
1110 nd = ir->opts.nrdf[i];
1111 reft = max(ir->opts.ref_t[i],0);
1116 if (IR_NVT_TROTTER(ir))
1118 /* contribution from the thermal momenta of the NH chain */
1122 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1123 /* contribution from the thermal variable of the NH chain */
1131 ener_npt += ndj*ixi[j]*kT;
1135 else /* Other non Trotter temperature NH control -- no chains yet. */
1137 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1138 ener_npt += nd*ixi[0]*kT;
1146 static real vrescale_gamdev(int ia, gmx_rng_t rng)
1147 /* Gamma distribution, adapted from numerical recipes */
1150 real am,e,s,v1,v2,x,y;
1157 for(j=1; j<=ia; j++)
1159 x *= gmx_rng_uniform_real(rng);
1173 v1 = gmx_rng_uniform_real(rng);
1174 v2 = 2.0*gmx_rng_uniform_real(rng)-1.0;
1176 while (v1*v1 + v2*v2 > 1.0 ||
1177 v1*v1*GMX_REAL_MAX < 3.0*ia);
1178 /* The last check above ensures that both x (3.0 > 2.0 in s)
1179 * and the pre-factor for e do not go out of range.
1183 s = sqrt(2.0*am + 1.0);
1187 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1189 while (gmx_rng_uniform_real(rng) > e);
1195 static real vrescale_sumnoises(int nn,gmx_rng_t rng)
1198 * Returns the sum of n independent gaussian noises squared
1199 * (i.e. equivalent to summing the square of the return values
1200 * of nn calls to gmx_rng_gaussian_real).xs
1206 } else if (nn == 1) {
1207 rr = gmx_rng_gaussian_real(rng);
1209 } else if (nn % 2 == 0) {
1210 return 2.0*vrescale_gamdev(nn/2,rng);
1212 rr = gmx_rng_gaussian_real(rng);
1213 return 2.0*vrescale_gamdev((nn-1)/2,rng) + rr*rr;
1217 static real vrescale_resamplekin(real kk,real sigma, int ndeg, real taut,
1221 * Generates a new value for the kinetic energy,
1222 * according to Bussi et al JCP (2007), Eq. (A7)
1223 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1224 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1225 * ndeg: number of degrees of freedom of the atoms to be thermalized
1226 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1231 factor = exp(-1.0/taut);
1235 rr = gmx_rng_gaussian_real(rng);
1238 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1,rng) + rr*rr)/ndeg - kk) +
1239 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1242 void vrescale_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt,
1243 double therm_integral[],gmx_rng_t rng)
1247 real Ek,Ek_ref1,Ek_ref,Ek_new;
1251 for(i=0; (i<opts->ngtc); i++)
1255 Ek = trace(ekind->tcstat[i].ekinf);
1259 Ek = trace(ekind->tcstat[i].ekinh);
1262 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1264 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1265 Ek_ref = Ek_ref1*opts->nrdf[i];
1267 Ek_new = vrescale_resamplekin(Ek,Ek_ref,opts->nrdf[i],
1268 opts->tau_t[i]/dt,rng);
1270 /* Analytically Ek_new>=0, but we check for rounding errors */
1273 ekind->tcstat[i].lambda = 0.0;
1277 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1280 therm_integral[i] -= Ek_new - Ek;
1284 fprintf(debug,"TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1285 i,Ek_ref,Ek,Ek_new,ekind->tcstat[i].lambda);
1290 ekind->tcstat[i].lambda = 1.0;
1295 real vrescale_energy(t_grpopts *opts,double therm_integral[])
1301 for(i=0; i<opts->ngtc; i++) {
1302 ener += therm_integral[i];
1308 void rescale_velocities(gmx_ekindata_t *ekind,t_mdatoms *mdatoms,
1309 int start,int end,rvec v[])
1312 t_grp_tcstat *tcstat;
1313 unsigned short *cACC,*cTC;
1318 tcstat = ekind->tcstat;
1323 gstat = ekind->grpstat;
1324 cACC = mdatoms->cACC;
1328 for(n=start; n<end; n++)
1338 /* Only scale the velocity component relative to the COM velocity */
1339 rvec_sub(v[n],gstat[ga].u,vrel);
1340 lg = tcstat[gt].lambda;
1341 for(d=0; d<DIM; d++)
1343 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1350 for(n=start; n<end; n++)
1356 lg = tcstat[gt].lambda;
1357 for(d=0; d<DIM; d++)
1366 /* set target temperatures if we are annealing */
1367 void update_annealing_target_temp(t_grpopts *opts,real t)
1370 real pert,thist=0,x;
1372 for(i=0;i<opts->ngtc;i++) {
1373 npoints = opts->anneal_npoints[i];
1374 switch (opts->annealing[i]) {
1378 /* calculate time modulo the period */
1379 pert = opts->anneal_time[i][npoints-1];
1381 thist = t - n*pert; /* modulo time */
1382 /* Make sure rounding didn't get us outside the interval */
1383 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1390 gmx_fatal(FARGS,"Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",i,opts->ngtc,npoints);
1392 /* We are doing annealing for this group if we got here,
1393 * and we have the (relative) time as thist.
1394 * calculate target temp */
1396 while ((j < npoints-1) && (thist>(opts->anneal_time[i][j+1])))
1398 if (j < npoints-1) {
1399 /* Found our position between points j and j+1.
1400 * Interpolate: x is the amount from j+1, (1-x) from point j
1401 * First treat possible jumps in temperature as a special case.
1403 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1404 opts->ref_t[i]=opts->anneal_temp[i][j+1];
1406 x = ((thist-opts->anneal_time[i][j])/
1407 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1408 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1412 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];