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) {
451 fprintf(fplog,"\nStep %s Warning: Pressure scaling more than 1%%.\n",
452 gmx_step_str(step,buf));
456 preserve_box_shape(ir,box_rel,boxv);
458 mtmul(boxv,box,t1); /* t1=boxv * b' */
462 /* Determine the scaling matrix mu for the coordinates */
465 t1[d][n] = box[d][n] + dt*boxv[d][n];
466 preserve_box_shape(ir,box_rel,t1);
467 /* t1 is the box at t+dt, determine mu as the relative change */
468 mmul_ur0(invbox,t1,mu);
471 void berendsen_pcoupl(FILE *fplog,gmx_large_int_t step,
472 t_inputrec *ir,real dt, tensor pres,matrix box,
476 real scalar_pressure, xy_pressure, p_corr_z;
477 char *ptr,buf[STRLEN];
480 * Calculate the scaling matrix mu
484 for(d=0; d<DIM; d++) {
485 scalar_pressure += pres[d][d]/DIM;
487 xy_pressure += pres[d][d]/(DIM-1);
489 /* Pressure is now in bar, everywhere. */
490 #define factor(d,m) (ir->compress[d][m]*dt/ir->tau_p)
492 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
493 * necessary for triclinic scaling
500 mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d] - scalar_pressure) /DIM;
503 case epctSEMIISOTROPIC:
505 mu[d][d] = 1.0 - factor(d,d)*(ir->ref_p[d][d]-xy_pressure)/DIM;
507 1.0 - factor(ZZ,ZZ)*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ])/DIM;
509 case epctANISOTROPIC:
512 mu[d][n] = (d==n ? 1.0 : 0.0)
513 -factor(d,n)*(ir->ref_p[d][n] - pres[d][n])/DIM;
515 case epctSURFACETENSION:
516 /* ir->ref_p[0/1] is the reference surface-tension times *
517 * the number of surfaces */
518 if (ir->compress[ZZ][ZZ])
519 p_corr_z = dt/ir->tau_p*(ir->ref_p[ZZ][ZZ] - pres[ZZ][ZZ]);
521 /* when the compressibity is zero, set the pressure correction *
522 * in the z-direction to zero to get the correct surface tension */
524 mu[ZZ][ZZ] = 1.0 - ir->compress[ZZ][ZZ]*p_corr_z;
525 for(d=0; d<DIM-1; d++)
526 mu[d][d] = 1.0 + factor(d,d)*(ir->ref_p[d][d]/(mu[ZZ][ZZ]*box[ZZ][ZZ])
527 - (pres[ZZ][ZZ]+p_corr_z - xy_pressure))/(DIM-1);
530 gmx_fatal(FARGS,"Berendsen pressure coupling type %s not supported yet\n",
531 EPCOUPLTYPETYPE(ir->epct));
534 /* To fullfill the orientation restrictions on triclinic boxes
535 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
536 * the other elements of mu to first order.
538 mu[YY][XX] += mu[XX][YY];
539 mu[ZZ][XX] += mu[XX][ZZ];
540 mu[ZZ][YY] += mu[YY][ZZ];
546 pr_rvecs(debug,0,"PC: pres ",pres,3);
547 pr_rvecs(debug,0,"PC: mu ",mu,3);
550 if (mu[XX][XX]<0.99 || mu[XX][XX]>1.01 ||
551 mu[YY][YY]<0.99 || mu[YY][YY]>1.01 ||
552 mu[ZZ][ZZ]<0.99 || mu[ZZ][ZZ]>1.01) {
554 sprintf(buf,"\nStep %s Warning: pressure scaling more than 1%%, "
556 gmx_step_str(step,buf2),mu[XX][XX],mu[YY][YY],mu[ZZ][ZZ]);
558 fprintf(fplog,"%s",buf);
559 fprintf(stderr,"%s",buf);
563 void berendsen_pscale(t_inputrec *ir,matrix mu,
564 matrix box,matrix box_rel,
565 int start,int nr_atoms,
566 rvec x[],unsigned short cFREEZE[],
569 ivec *nFreeze=ir->opts.nFreeze;
572 /* Scale the positions */
573 for (n=start; n<start+nr_atoms; n++) {
578 x[n][XX] = mu[XX][XX]*x[n][XX]+mu[YY][XX]*x[n][YY]+mu[ZZ][XX]*x[n][ZZ];
580 x[n][YY] = mu[YY][YY]*x[n][YY]+mu[ZZ][YY]*x[n][ZZ];
582 x[n][ZZ] = mu[ZZ][ZZ]*x[n][ZZ];
584 /* compute final boxlengths */
585 for (d=0; d<DIM; d++) {
586 box[d][XX] = mu[XX][XX]*box[d][XX]+mu[YY][XX]*box[d][YY]+mu[ZZ][XX]*box[d][ZZ];
587 box[d][YY] = mu[YY][YY]*box[d][YY]+mu[ZZ][YY]*box[d][ZZ];
588 box[d][ZZ] = mu[ZZ][ZZ]*box[d][ZZ];
591 preserve_box_shape(ir,box_rel,box);
593 /* (un)shifting should NOT be done after this,
594 * since the box vectors might have changed
596 inc_nrnb(nrnb,eNR_PCOUPL,nr_atoms);
599 void berendsen_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt)
607 for(i=0; (i<opts->ngtc); i++)
611 T = ekind->tcstat[i].T;
615 T = ekind->tcstat[i].Th;
618 if ((opts->tau_t[i] > 0) && (T > 0.0)) {
620 reft = max(0.0,opts->ref_t[i]);
621 lll = sqrt(1.0 + (dt/opts->tau_t[i])*(reft/T-1.0));
622 ekind->tcstat[i].lambda = max(min(lll,1.25),0.8);
625 ekind->tcstat[i].lambda = 1.0;
629 fprintf(debug,"TC: group %d: T: %g, Lambda: %g\n",
630 i,T,ekind->tcstat[i].lambda);
634 void nosehoover_tcoupl(t_grpopts *opts,gmx_ekindata_t *ekind,real dt,
635 double xi[],double vxi[], t_extmass *MassQ)
640 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
642 for(i=0; (i<opts->ngtc); i++) {
643 reft = max(0.0,opts->ref_t[i]);
645 vxi[i] += dt*MassQ->Qinv[i]*(ekind->tcstat[i].Th - reft);
646 xi[i] += dt*(oldvxi + vxi[i])*0.5;
650 t_state *init_bufstate(const t_state *template_state)
653 int nc = template_state->nhchainlength;
655 snew(state->nosehoover_xi,nc*template_state->ngtc);
656 snew(state->nosehoover_vxi,nc*template_state->ngtc);
657 snew(state->therm_integral,template_state->ngtc);
658 snew(state->nhpres_xi,nc*template_state->nnhpres);
659 snew(state->nhpres_vxi,nc*template_state->nnhpres);
664 void destroy_bufstate(t_state *state)
668 sfree(state->nosehoover_xi);
669 sfree(state->nosehoover_vxi);
670 sfree(state->therm_integral);
671 sfree(state->nhpres_xi);
672 sfree(state->nhpres_vxi);
676 void trotter_update(t_inputrec *ir,gmx_large_int_t step, gmx_ekindata_t *ekind,
677 gmx_enerdata_t *enerd, t_state *state,
678 tensor vir, t_mdatoms *md,
679 t_extmass *MassQ, int **trotter_seqlist, int trotter_seqno)
682 int n,i,j,d,ntgrp,ngtc,gc=0;
683 t_grp_tcstat *tcstat;
685 gmx_large_int_t step_eff;
686 real ecorr,pcorr,dvdlcorr;
687 real bmass,qmass,reft,kT,dt,nd;
688 tensor dumpres,dumvir;
689 double *scalefac,dtc;
694 if (trotter_seqno <= ettTSEQ2)
696 step_eff = step-1; /* the velocity verlet calls are actually out of order -- the first half step
697 is actually the last half step from the previous step. Thus the first half step
698 actually corresponds to the n-1 step*/
704 bCouple = (ir->nsttcouple == 1 ||
705 do_per_step(step_eff+ir->nsttcouple,ir->nsttcouple));
707 trotter_seq = trotter_seqlist[trotter_seqno];
709 /* signal we are returning if nothing is going to be done in this routine */
710 if ((trotter_seq[0] == etrtSKIPALL) || !(bCouple))
715 dtc = ir->nsttcouple*ir->delta_t;
716 opts = &(ir->opts); /* just for ease of referencing */
718 snew(scalefac,opts->ngtc);
723 /* execute the series of trotter updates specified in the trotterpart array */
725 for (i=0;i<NTROTTERPARTS;i++){
726 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
727 if ((trotter_seq[i] == etrtBAROV2) || (trotter_seq[i] == etrtBARONHC2) || (trotter_seq[i] == etrtNHC2))
736 switch (trotter_seq[i])
740 boxv_trotter(ir,&(state->veta),dt,state->box,ekind,vir,
741 enerd->term[F_PDISPCORR],enerd->term[F_DISPCORR],MassQ);
745 NHC_trotter(opts,state->nnhpres,ekind,dt,state->nhpres_xi,
746 state->nhpres_vxi,NULL,&(state->veta),MassQ,FALSE);
750 NHC_trotter(opts,opts->ngtc,ekind,dt,state->nosehoover_xi,
751 state->nosehoover_vxi,scalefac,NULL,MassQ,(ir->eI==eiVV));
752 /* need to rescale the kinetic energies and velocities here. Could
753 scale the velocities later, but we need them scaled in order to
754 produce the correct outputs, so we'll scale them here. */
756 for (i=0; i<ngtc;i++)
758 tcstat = &ekind->tcstat[i];
759 tcstat->vscale_nhc = scalefac[i];
760 tcstat->ekinscaleh_nhc *= (scalefac[i]*scalefac[i]);
761 tcstat->ekinscalef_nhc *= (scalefac[i]*scalefac[i]);
763 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
764 /* but do we actually need the total? */
766 /* modify the velocities as well */
767 for (n=md->start;n<md->start+md->homenr;n++)
775 state->v[n][d] *= scalefac[gc];
782 sumv[d] += (state->v[n][d])/md->invmass[n];
791 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
799 consk[d] = sumv[d]*exp((1 + 1.0/opts->nrdf[0])*((1.0/DIM)*log(det(state->box)/state->vol0)) + state->nosehoover_xi[0]);
801 fprintf(debug,"Conserved kappa: %15.8f %15.8f %15.8f\n",consk[0],consk[1],consk[2]);
808 int **init_npt_vars(t_inputrec *ir, t_state *state, t_extmass *MassQ, gmx_bool bTrotter)
810 int n,i,j,d,ntgrp,ngtc,nnhpres,nh,gc=0;
811 t_grp_tcstat *tcstat;
813 real ecorr,pcorr,dvdlcorr;
814 real bmass,qmass,reft,kT,dt,ndj,nd;
815 tensor dumpres,dumvir;
818 opts = &(ir->opts); /* just for ease of referencing */
820 nnhpres = state->nnhpres;
821 nh = state->nhchainlength;
823 if (ir->eI == eiMD) {
824 snew(MassQ->Qinv,ngtc);
825 for(i=0; (i<ngtc); i++)
827 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
829 MassQ->Qinv[i]=1.0/(sqr(opts->tau_t[i]/M_2PI)*opts->ref_t[i]);
837 else if (EI_VV(ir->eI))
839 /* Set pressure variables */
841 if (state->vol0 == 0)
843 state->vol0 = det(state->box); /* because we start by defining a fixed compressibility,
844 we need the volume at this compressibility to solve the problem */
847 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
848 /* Investigate this more -- is this the right mass to make it? */
849 MassQ->Winv = (PRESFAC*trace(ir->compress)*BOLTZ*opts->ref_t[0])/(DIM*state->vol0*sqr(ir->tau_p/M_2PI));
850 /* An alternate mass definition, from Tuckerman et al. */
851 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
856 MassQ->Winvm[d][n]= PRESFAC*ir->compress[d][n]/(state->vol0*sqr(ir->tau_p/M_2PI));
857 /* not clear this is correct yet for the anisotropic case*/
860 /* Allocate space for thermostat variables */
861 snew(MassQ->Qinv,ngtc*nh);
863 /* now, set temperature variables */
864 for(i=0; i<ngtc; i++)
866 if ((opts->tau_t[i] > 0) && (opts->ref_t[i] > 0))
868 reft = max(0.0,opts->ref_t[i]);
881 MassQ->Qinv[i*nh+j] = 1.0/(sqr(opts->tau_t[i]/M_2PI)*ndj*kT);
889 MassQ->Qinv[i*nh+j] = 0.0;
895 /* first, initialize clear all the trotter calls */
896 snew(trotter_seq,ettTSEQMAX);
897 for (i=0;i<ettTSEQMAX;i++)
899 snew(trotter_seq[i],NTROTTERPARTS);
900 for (j=0;j<NTROTTERPARTS;j++) {
901 trotter_seq[i][j] = etrtNONE;
903 trotter_seq[i][0] = etrtSKIPALL;
908 /* no trotter calls, so we never use the values in the array.
909 * We access them (so we need to define them, but ignore
915 /* compute the kinetic energy by using the half step velocities or
916 * the kinetic energies, depending on the order of the trotter calls */
920 if (IR_NPT_TROTTER(ir))
922 /* This is the complicated version - there are 4 possible calls, depending on ordering.
923 We start with the initial one. */
924 /* first, a round that estimates veta. */
925 trotter_seq[0][0] = etrtBAROV;
927 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
929 /* The first half trotter update */
930 trotter_seq[2][0] = etrtBAROV;
931 trotter_seq[2][1] = etrtNHC;
932 trotter_seq[2][2] = etrtBARONHC;
934 /* The second half trotter update */
935 trotter_seq[3][0] = etrtBARONHC;
936 trotter_seq[3][1] = etrtNHC;
937 trotter_seq[3][2] = etrtBAROV;
939 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
944 if (IR_NVT_TROTTER(ir))
946 /* This is the easy version - there are only two calls, both the same.
947 Otherwise, even easier -- no calls */
948 trotter_seq[2][0] = etrtNHC;
949 trotter_seq[3][0] = etrtNHC;
952 } else if (ir->eI==eiVVAK) {
953 if (IR_NPT_TROTTER(ir))
955 /* This is the complicated version - there are 4 possible calls, depending on ordering.
956 We start with the initial one. */
957 /* first, a round that estimates veta. */
958 trotter_seq[0][0] = etrtBAROV;
960 /* The first half trotter update, part 1 -- double update, because it commutes */
961 trotter_seq[1][0] = etrtNHC;
963 /* The first half trotter update, part 2 */
964 trotter_seq[2][0] = etrtBAROV;
965 trotter_seq[2][1] = etrtBARONHC;
967 /* The second half trotter update, part 1 */
968 trotter_seq[3][0] = etrtBARONHC;
969 trotter_seq[3][1] = etrtBAROV;
971 /* The second half trotter update -- blank for now */
972 trotter_seq[4][0] = etrtNHC;
976 if (IR_NVT_TROTTER(ir))
978 /* This is the easy version - there is only one call, both the same.
979 Otherwise, even easier -- no calls */
980 trotter_seq[1][0] = etrtNHC;
981 trotter_seq[4][0] = etrtNHC;
990 bmass = DIM*DIM; /* recommended mass parameters for isotropic barostat */
993 snew(MassQ->QPinv,nnhpres*opts->nhchainlength);
995 /* barostat temperature */
996 if ((ir->tau_p > 0) && (opts->ref_t[0] > 0))
998 reft = max(0.0,opts->ref_t[0]);
1000 for (i=0;i<nnhpres;i++) {
1010 MassQ->QPinv[i*opts->nhchainlength+j] = 1.0/(sqr(opts->tau_t[0]/M_2PI)*qmass*kT);
1016 for (i=0;i<nnhpres;i++) {
1019 MassQ->QPinv[i*nh+j]=0.0;
1026 real NPT_energy(t_inputrec *ir, t_state *state, t_extmass *MassQ)
1028 int i,j,nd,ndj,bmass,qmass,ngtcall;
1029 real ener_npt,reft,eta,kT,tau;
1033 int nh = state->nhchainlength;
1037 /* now we compute the contribution of the pressure to the conserved quantity*/
1039 if (ir->epc==epcMTTK)
1041 /* find the volume, and the kinetic energy of the volume */
1046 /* contribution from the pressure momenenta */
1047 ener_npt += 0.5*sqr(state->veta)/MassQ->Winv;
1049 /* contribution from the PV term */
1050 vol = det(state->box);
1051 ener_npt += vol*trace(ir->ref_p)/(DIM*PRESFAC);
1054 case epctANISOTROPIC:
1058 case epctSURFACETENSION:
1061 case epctSEMIISOTROPIC:
1069 if (IR_NPT_TROTTER(ir))
1071 /* add the energy from the barostat thermostat chain */
1072 for (i=0;i<state->nnhpres;i++) {
1074 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1075 ivxi = &state->nhpres_vxi[i*nh];
1076 ixi = &state->nhpres_xi[i*nh];
1077 iQinv = &(MassQ->QPinv[i*nh]);
1078 reft = max(ir->opts.ref_t[0],0); /* using 'System' temperature */
1085 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1086 /* contribution from the thermal variable of the NH chain */
1087 ener_npt += ixi[j]*kT;
1091 fprintf(debug,"P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f",i,j,ivxi[j],ixi[j]);
1099 for(i=0; i<ir->opts.ngtc; i++)
1101 ixi = &state->nosehoover_xi[i*nh];
1102 ivxi = &state->nosehoover_vxi[i*nh];
1103 iQinv = &(MassQ->Qinv[i*nh]);
1105 nd = ir->opts.nrdf[i];
1106 reft = max(ir->opts.ref_t[i],0);
1111 if (IR_NVT_TROTTER(ir))
1113 /* contribution from the thermal momenta of the NH chain */
1117 ener_npt += 0.5*sqr(ivxi[j])/iQinv[j];
1118 /* contribution from the thermal variable of the NH chain */
1126 ener_npt += ndj*ixi[j]*kT;
1130 else /* Other non Trotter temperature NH control -- no chains yet. */
1132 ener_npt += 0.5*BOLTZ*nd*sqr(ivxi[0])/iQinv[0];
1133 ener_npt += nd*ixi[0]*kT;
1141 static real vrescale_gamdev(int ia, gmx_rng_t rng)
1142 /* Gamma distribution, adapted from numerical recipes */
1145 real am,e,s,v1,v2,x,y;
1152 for(j=1; j<=ia; j++)
1154 x *= gmx_rng_uniform_real(rng);
1168 v1 = gmx_rng_uniform_real(rng);
1169 v2 = 2.0*gmx_rng_uniform_real(rng)-1.0;
1171 while (v1*v1 + v2*v2 > 1.0 ||
1172 v1*v1*GMX_REAL_MAX < 3.0*ia);
1173 /* The last check above ensures that both x (3.0 > 2.0 in s)
1174 * and the pre-factor for e do not go out of range.
1178 s = sqrt(2.0*am + 1.0);
1182 e = (1.0 + y*y)*exp(am*log(x/am) - s*y);
1184 while (gmx_rng_uniform_real(rng) > e);
1190 static real vrescale_sumnoises(int nn,gmx_rng_t rng)
1193 * Returns the sum of n independent gaussian noises squared
1194 * (i.e. equivalent to summing the square of the return values
1195 * of nn calls to gmx_rng_gaussian_real).xs
1201 } else if (nn == 1) {
1202 rr = gmx_rng_gaussian_real(rng);
1204 } else if (nn % 2 == 0) {
1205 return 2.0*vrescale_gamdev(nn/2,rng);
1207 rr = gmx_rng_gaussian_real(rng);
1208 return 2.0*vrescale_gamdev((nn-1)/2,rng) + rr*rr;
1212 static real vrescale_resamplekin(real kk,real sigma, int ndeg, real taut,
1216 * Generates a new value for the kinetic energy,
1217 * according to Bussi et al JCP (2007), Eq. (A7)
1218 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1219 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1220 * ndeg: number of degrees of freedom of the atoms to be thermalized
1221 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1226 factor = exp(-1.0/taut);
1230 rr = gmx_rng_gaussian_real(rng);
1233 (1.0 - factor)*(sigma*(vrescale_sumnoises(ndeg-1,rng) + rr*rr)/ndeg - kk) +
1234 2.0*rr*sqrt(kk*sigma/ndeg*(1.0 - factor)*factor);
1237 void vrescale_tcoupl(t_inputrec *ir,gmx_ekindata_t *ekind,real dt,
1238 double therm_integral[],gmx_rng_t rng)
1242 real Ek,Ek_ref1,Ek_ref,Ek_new;
1246 for(i=0; (i<opts->ngtc); i++)
1250 Ek = trace(ekind->tcstat[i].ekinf);
1254 Ek = trace(ekind->tcstat[i].ekinh);
1257 if (opts->tau_t[i] >= 0 && opts->nrdf[i] > 0 && Ek > 0)
1259 Ek_ref1 = 0.5*opts->ref_t[i]*BOLTZ;
1260 Ek_ref = Ek_ref1*opts->nrdf[i];
1262 Ek_new = vrescale_resamplekin(Ek,Ek_ref,opts->nrdf[i],
1263 opts->tau_t[i]/dt,rng);
1265 /* Analytically Ek_new>=0, but we check for rounding errors */
1268 ekind->tcstat[i].lambda = 0.0;
1272 ekind->tcstat[i].lambda = sqrt(Ek_new/Ek);
1275 therm_integral[i] -= Ek_new - Ek;
1279 fprintf(debug,"TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1280 i,Ek_ref,Ek,Ek_new,ekind->tcstat[i].lambda);
1285 ekind->tcstat[i].lambda = 1.0;
1290 real vrescale_energy(t_grpopts *opts,double therm_integral[])
1296 for(i=0; i<opts->ngtc; i++) {
1297 ener += therm_integral[i];
1303 void rescale_velocities(gmx_ekindata_t *ekind,t_mdatoms *mdatoms,
1304 int start,int end,rvec v[])
1307 t_grp_tcstat *tcstat;
1308 unsigned short *cACC,*cTC;
1313 tcstat = ekind->tcstat;
1318 gstat = ekind->grpstat;
1319 cACC = mdatoms->cACC;
1323 for(n=start; n<end; n++)
1333 /* Only scale the velocity component relative to the COM velocity */
1334 rvec_sub(v[n],gstat[ga].u,vrel);
1335 lg = tcstat[gt].lambda;
1336 for(d=0; d<DIM; d++)
1338 v[n][d] = gstat[ga].u[d] + lg*vrel[d];
1345 for(n=start; n<end; n++)
1351 lg = tcstat[gt].lambda;
1352 for(d=0; d<DIM; d++)
1361 /* set target temperatures if we are annealing */
1362 void update_annealing_target_temp(t_grpopts *opts,real t)
1365 real pert,thist=0,x;
1367 for(i=0;i<opts->ngtc;i++) {
1368 npoints = opts->anneal_npoints[i];
1369 switch (opts->annealing[i]) {
1373 /* calculate time modulo the period */
1374 pert = opts->anneal_time[i][npoints-1];
1376 thist = t - n*pert; /* modulo time */
1377 /* Make sure rounding didn't get us outside the interval */
1378 if (fabs(thist-pert) < GMX_REAL_EPS*100)
1385 gmx_fatal(FARGS,"Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",i,opts->ngtc,npoints);
1387 /* We are doing annealing for this group if we got here,
1388 * and we have the (relative) time as thist.
1389 * calculate target temp */
1391 while ((j < npoints-1) && (thist>(opts->anneal_time[i][j+1])))
1393 if (j < npoints-1) {
1394 /* Found our position between points j and j+1.
1395 * Interpolate: x is the amount from j+1, (1-x) from point j
1396 * First treat possible jumps in temperature as a special case.
1398 if ((opts->anneal_time[i][j+1]-opts->anneal_time[i][j]) < GMX_REAL_EPS*100)
1399 opts->ref_t[i]=opts->anneal_temp[i][j+1];
1401 x = ((thist-opts->anneal_time[i][j])/
1402 (opts->anneal_time[i][j+1]-opts->anneal_time[i][j]));
1403 opts->ref_t[i] = x*opts->anneal_temp[i][j+1]+(1-x)*opts->anneal_temp[i][j];
1407 opts->ref_t[i] = opts->anneal_temp[i][npoints-1];