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 NTROTTERCALLS 5
54 #define NTROTTERPARTS 3
56 /* these integration routines are only referenced inside this file */
57 static void NHC_trotter(t_grpopts
*opts
,int nvar
, gmx_ekindata_t
*ekind
,real dtfull
,
58 double xi
[],double vxi
[], double scalefac
[], real
*veta
, t_extmass
*MassQ
, bool bEkinAveVel
)
61 /* general routine for both barostat and thermostat nose hoover chains */
63 int i
,j
,mi
,mj
,jmax
,nd
;
64 double Ekin
,Efac
,reft
,kT
;
72 int ns
= SUZUKI_YOSHIDA_NUM
; /* set the degree of integration in the types/state.h file */
73 int nh
= opts
->nhchainlength
;
76 mstepsi
= mstepsj
= ns
;
78 /* if scalefac is NULL, we are doing the NHC of the barostat */
81 if (scalefac
== NULL
) {
85 for (i
=0; i
<nvar
; i
++)
88 /* make it easier to iterate by selecting
89 out the sub-array that corresponds to this T group */
94 iQinv
= &(MassQ
->QPinv
[i
*nh
]);
95 nd
= 1; /* THIS WILL CHANGE IF NOT ISOTROPIC */
96 reft
= max(0.0,opts
->ref_t
[0]);
97 Ekin
= sqr(*veta
)/MassQ
->Winv
;
99 iQinv
= &(MassQ
->Qinv
[i
*nh
]);
100 tcstat
= &ekind
->tcstat
[i
];
102 reft
= max(0.0,opts
->ref_t
[i
]);
105 Ekin
= 2*trace(tcstat
->ekinf
)*tcstat
->ekinscalef_nhc
;
107 Ekin
= 2*trace(tcstat
->ekinh
)*tcstat
->ekinscaleh_nhc
;
112 for(mi
=0;mi
<mstepsi
;mi
++)
114 for(mj
=0;mj
<mstepsj
;mj
++)
116 /* weighting for this step using Suzuki-Yoshida integration - fixed at 5 */
117 dt
= sy_const
[ns
][mj
] * dtfull
/ mstepsi
;
119 /* compute the thermal forces */
120 GQ
[0] = iQinv
[0]*(Ekin
- nd
*kT
);
124 if (iQinv
[j
+1] > 0) {
125 /* we actually don't need to update here if we save the
126 state of the GQ, but it's easier to just recompute*/
127 GQ
[j
+1] = iQinv
[j
+1]*((sqr(ivxi
[j
])/iQinv
[j
])-kT
);
133 ivxi
[nh
-1] += 0.25*dt
*GQ
[nh
-1];
136 Efac
= exp(-0.125*dt
*ivxi
[j
]);
137 ivxi
[j
-1] = Efac
*(ivxi
[j
-1]*Efac
+ 0.25*dt
*GQ
[j
-1]);
140 Efac
= exp(-0.5*dt
*ivxi
[0]);
148 /* Issue - if the KE is an average of the last and the current temperatures, then we might not be
149 able to scale the kinetic energy directly with this factor. Might take more bookkeeping -- have to
150 think about this a bit more . . . */
152 GQ
[0] = iQinv
[0]*(Ekin
- nd
*kT
);
154 /* update thermostat positions */
157 ixi
[j
] += 0.5*dt
*ivxi
[j
];
162 Efac
= exp(-0.125*dt
*ivxi
[j
+1]);
163 ivxi
[j
] = Efac
*(ivxi
[j
]*Efac
+ 0.25*dt
*GQ
[j
]);
164 if (iQinv
[j
+1] > 0) {
165 GQ
[j
+1] = iQinv
[j
+1]*((sqr(ivxi
[j
])/iQinv
[j
])-kT
);
170 ivxi
[nh
-1] += 0.25*dt
*GQ
[nh
-1];
177 static void boxv_trotter(t_inputrec
*ir
, real
*veta
, real dt
, tensor box
,
178 gmx_ekindata_t
*ekind
, tensor vir
, real pcorr
, real ecorr
, t_extmass
*MassQ
)
185 tensor Winvm
,ekinmod
,localpres
;
187 /* The heat bath is coupled to a separate barostat, the last temperature group. In the
188 2006 Tuckerman et al paper., the order is iL_{T_baro} iL {T_part}
191 if (ir
->epct
==epctSEMIISOTROPIC
)
200 /* eta is in pure units. veta is in units of ps^-1. GW is in
201 units of ps^-2. However, eta has a reference of 1 nm^3, so care must be
202 taken to use only RATIOS of eta in updating the volume. */
204 /* we take the partial pressure tensors, modify the
205 kinetic energy tensor, and recovert to pressure */
207 if (ir
->opts
.nrdf
[0]==0)
209 gmx_fatal(FARGS
,"Barostat is coupled to a T-group with no degrees of freedom\n");
211 /* alpha factor for phase space volume, then multiply by the ekin scaling factor. */
212 alpha
= 1.0 + DIM
/((double)ir
->opts
.nrdf
[0]);
213 alpha
*= ekind
->tcstat
[0].ekinscalef_nhc
;
214 msmul(ekind
->ekin
,alpha
,ekinmod
);
216 /* for now, we use Elr = 0, because if you want to get it right, you
217 really should be using PME. Maybe print a warning? */
219 pscal
= calc_pres(ir
->ePBC
,nwall
,box
,ekinmod
,vir
,localpres
,0.0) + pcorr
;
222 GW
= (vol
*(MassQ
->Winv
/PRESFAC
))*(DIM
*pscal
- trace(ir
->ref_p
)); /* W is in ps^2 * bar * nm^3 */
228 * This file implements temperature and pressure coupling algorithms:
229 * For now only the Weak coupling and the modified weak coupling.
231 * Furthermore computation of pressure and temperature is done here
235 real
calc_pres(int ePBC
,int nwall
,matrix box
,tensor ekin
,tensor vir
,
236 tensor pres
,real Elr
)
241 if (ePBC
==epbcNONE
|| (ePBC
==epbcXY
&& nwall
!=2))
244 /* Uitzoeken welke ekin hier van toepassing is, zie Evans & Morris - E.
245 * Wrs. moet de druktensor gecorrigeerd worden voor de netto stroom in
249 /* Long range correction for periodic systems, see
251 * divide by 6 because it is multiplied by fac later on.
252 * If Elr = 0, no correction is made.
255 /* This formula should not be used with Ewald or PME,
256 * where the full long-range virial is calculated. EL 990823
260 fac
=PRESFAC
*2.0/det(box
);
261 for(n
=0; (n
<DIM
); n
++)
262 for(m
=0; (m
<DIM
); m
++)
263 pres
[n
][m
]=(ekin
[n
][m
]-vir
[n
][m
]+Plr
)*fac
;
266 pr_rvecs(debug
,0,"PC: pres",pres
,DIM
);
267 pr_rvecs(debug
,0,"PC: ekin",ekin
,DIM
);
268 pr_rvecs(debug
,0,"PC: vir ",vir
, DIM
);
269 pr_rvecs(debug
,0,"PC: box ",box
, DIM
);
272 return trace(pres
)/DIM
;
275 real
calc_temp(real ekin
,real nrdf
)
278 return (2.0*ekin
)/(nrdf
*BOLTZ
);
283 void parrinellorahman_pcoupl(FILE *fplog
,gmx_large_int_t step
,
284 t_inputrec
*ir
,real dt
,tensor pres
,
285 tensor box
,tensor box_rel
,tensor boxv
,
286 tensor M
,matrix mu
,bool bFirstStep
)
288 /* This doesn't do any coordinate updating. It just
289 * integrates the box vector equations from the calculated
290 * acceleration due to pressure difference. We also compute
291 * the tensor M which is used in update to couple the particle
292 * coordinates to the box vectors.
294 * In Nose and Klein (Mol.Phys 50 (1983) no 5., p 1055) this is
297 * M_nk = (h') * (h' * h + h' h) * h
299 * with the dots denoting time derivatives and h is the transformation from
300 * the scaled frame to the real frame, i.e. the TRANSPOSE of the box.
301 * This also goes for the pressure and M tensors - they are transposed relative
302 * to ours. Our equation thus becomes:
305 * M_gmx = M_nk' = b * (b * b' + b * b') * b'
307 * where b is the gromacs box matrix.
308 * Our box accelerations are given by
310 * b = vol/W inv(box') * (P-ref_P) (=h')
315 real vol
=box
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
316 real atot
,arel
,change
,maxchange
,xy_pressure
;
317 tensor invbox
,pdiff
,t1
,t2
;
321 m_inv_ur0(box
,invbox
);
324 /* Note that PRESFAC does not occur here.
325 * The pressure and compressibility always occur as a product,
326 * therefore the pressure unit drops out.
328 maxl
=max(box
[XX
][XX
],box
[YY
][YY
]);
329 maxl
=max(maxl
,box
[ZZ
][ZZ
]);
333 (4*M_PI
*M_PI
*ir
->compress
[d
][n
])/(3*ir
->tau_p
*ir
->tau_p
*maxl
);
335 m_sub(pres
,ir
->ref_p
,pdiff
);
337 if(ir
->epct
==epctSURFACETENSION
) {
338 /* Unlike Berendsen coupling it might not be trivial to include a z
339 * pressure correction here? On the other hand we don't scale the
340 * box momentarily, but change accelerations, so it might not be crucial.
342 xy_pressure
=0.5*(pres
[XX
][XX
]+pres
[YY
][YY
]);
344 pdiff
[d
][d
]=(xy_pressure
-(pres
[ZZ
][ZZ
]-ir
->ref_p
[d
][d
]/box
[d
][d
]));
347 tmmul(invbox
,pdiff
,t1
);
348 /* Move the off-diagonal elements of the 'force' to one side to ensure
349 * that we obey the box constraints.
353 t1
[d
][n
] += t1
[n
][d
];
359 case epctANISOTROPIC
:
362 t1
[d
][n
] *= winv
[d
][n
]*vol
;
365 /* calculate total volume acceleration */
366 atot
=box
[XX
][XX
]*box
[YY
][YY
]*t1
[ZZ
][ZZ
]+
367 box
[XX
][XX
]*t1
[YY
][YY
]*box
[ZZ
][ZZ
]+
368 t1
[XX
][XX
]*box
[YY
][YY
]*box
[ZZ
][ZZ
];
370 /* set all RELATIVE box accelerations equal, and maintain total V
374 t1
[d
][n
] = winv
[0][0]*vol
*arel
*box
[d
][n
];
376 case epctSEMIISOTROPIC
:
377 case epctSURFACETENSION
:
378 /* Note the correction to pdiff above for surftens. coupling */
380 /* calculate total XY volume acceleration */
381 atot
=box
[XX
][XX
]*t1
[YY
][YY
]+t1
[XX
][XX
]*box
[YY
][YY
];
382 arel
=atot
/(2*box
[XX
][XX
]*box
[YY
][YY
]);
383 /* set RELATIVE XY box accelerations equal, and maintain total V
384 * change speed. Dont change the third box vector accelerations */
387 t1
[d
][n
] = winv
[d
][n
]*vol
*arel
*box
[d
][n
];
389 t1
[ZZ
][n
] *= winv
[d
][n
]*vol
;
392 gmx_fatal(FARGS
,"Parrinello-Rahman pressure coupling type %s "
393 "not supported yet\n",EPCOUPLTYPETYPE(ir
->epct
));
400 boxv
[d
][n
] += dt
*t1
[d
][n
];
402 /* We do NOT update the box vectors themselves here, since
403 * we need them for shifting later. It is instead done last
404 * in the update() routine.
407 /* Calculate the change relative to diagonal elements-
408 since it's perfectly ok for the off-diagonal ones to
409 be zero it doesn't make sense to check the change relative
413 change
=fabs(dt
*boxv
[d
][n
]/box
[d
][d
]);
415 if (change
>maxchange
)
419 if (maxchange
> 0.01 && fplog
) {
421 fprintf(fplog
,"\nStep %s Warning: Pressure scaling more than 1%%.\n",
422 gmx_step_str(step
,buf
));
426 preserve_box_shape(ir
,box_rel
,boxv
);
428 mtmul(boxv
,box
,t1
); /* t1=boxv * b' */
432 /* Determine the scaling matrix mu for the coordinates */
435 t1
[d
][n
] = box
[d
][n
] + dt
*boxv
[d
][n
];
436 preserve_box_shape(ir
,box_rel
,t1
);
437 /* t1 is the box at t+dt, determine mu as the relative change */
438 mmul_ur0(invbox
,t1
,mu
);
441 void berendsen_pcoupl(FILE *fplog
,gmx_large_int_t step
,
442 t_inputrec
*ir
,real dt
, tensor pres
,matrix box
,
446 real scalar_pressure
, xy_pressure
, p_corr_z
;
447 char *ptr
,buf
[STRLEN
];
450 * Calculate the scaling matrix mu
454 for(d
=0; d
<DIM
; d
++) {
455 scalar_pressure
+= pres
[d
][d
]/DIM
;
457 xy_pressure
+= pres
[d
][d
]/(DIM
-1);
459 /* Pressure is now in bar, everywhere. */
460 #define factor(d,m) (ir->compress[d][m]*dt/ir->tau_p)
462 /* mu has been changed from pow(1+...,1/3) to 1+.../3, since this is
463 * necessary for triclinic scaling
470 mu
[d
][d
] = 1.0 - factor(d
,d
)*(ir
->ref_p
[d
][d
] - scalar_pressure
) /DIM
;
473 case epctSEMIISOTROPIC
:
475 mu
[d
][d
] = 1.0 - factor(d
,d
)*(ir
->ref_p
[d
][d
]-xy_pressure
)/DIM
;
477 1.0 - factor(ZZ
,ZZ
)*(ir
->ref_p
[ZZ
][ZZ
] - pres
[ZZ
][ZZ
])/DIM
;
479 case epctANISOTROPIC
:
482 mu
[d
][n
] = (d
==n
? 1.0 : 0.0)
483 -factor(d
,n
)*(ir
->ref_p
[d
][n
] - pres
[d
][n
])/DIM
;
485 case epctSURFACETENSION
:
486 /* ir->ref_p[0/1] is the reference surface-tension times *
487 * the number of surfaces */
488 if (ir
->compress
[ZZ
][ZZ
])
489 p_corr_z
= dt
/ir
->tau_p
*(ir
->ref_p
[ZZ
][ZZ
] - pres
[ZZ
][ZZ
]);
491 /* when the compressibity is zero, set the pressure correction *
492 * in the z-direction to zero to get the correct surface tension */
494 mu
[ZZ
][ZZ
] = 1.0 - ir
->compress
[ZZ
][ZZ
]*p_corr_z
;
495 for(d
=0; d
<DIM
-1; d
++)
496 mu
[d
][d
] = 1.0 + factor(d
,d
)*(ir
->ref_p
[d
][d
]/(mu
[ZZ
][ZZ
]*box
[ZZ
][ZZ
])
497 - (pres
[ZZ
][ZZ
]+p_corr_z
- xy_pressure
))/(DIM
-1);
500 gmx_fatal(FARGS
,"Berendsen pressure coupling type %s not supported yet\n",
501 EPCOUPLTYPETYPE(ir
->epct
));
504 /* To fullfill the orientation restrictions on triclinic boxes
505 * we will set mu_yx, mu_zx and mu_zy to 0 and correct
506 * the other elements of mu to first order.
508 mu
[YY
][XX
] += mu
[XX
][YY
];
509 mu
[ZZ
][XX
] += mu
[XX
][ZZ
];
510 mu
[ZZ
][YY
] += mu
[YY
][ZZ
];
516 pr_rvecs(debug
,0,"PC: pres ",pres
,3);
517 pr_rvecs(debug
,0,"PC: mu ",mu
,3);
520 if (mu
[XX
][XX
]<0.99 || mu
[XX
][XX
]>1.01 ||
521 mu
[YY
][YY
]<0.99 || mu
[YY
][YY
]>1.01 ||
522 mu
[ZZ
][ZZ
]<0.99 || mu
[ZZ
][ZZ
]>1.01) {
524 sprintf(buf
,"\nStep %s Warning: pressure scaling more than 1%%, "
526 gmx_step_str(step
,buf2
),mu
[XX
][XX
],mu
[YY
][YY
],mu
[ZZ
][ZZ
]);
528 fprintf(fplog
,"%s",buf
);
529 fprintf(stderr
,"%s",buf
);
533 void berendsen_pscale(t_inputrec
*ir
,matrix mu
,
534 matrix box
,matrix box_rel
,
535 int start
,int nr_atoms
,
536 rvec x
[],unsigned short cFREEZE
[],
539 ivec
*nFreeze
=ir
->opts
.nFreeze
;
542 /* Scale the positions */
543 for (n
=start
; n
<start
+nr_atoms
; n
++) {
548 x
[n
][XX
] = mu
[XX
][XX
]*x
[n
][XX
]+mu
[YY
][XX
]*x
[n
][YY
]+mu
[ZZ
][XX
]*x
[n
][ZZ
];
550 x
[n
][YY
] = mu
[YY
][YY
]*x
[n
][YY
]+mu
[ZZ
][YY
]*x
[n
][ZZ
];
552 x
[n
][ZZ
] = mu
[ZZ
][ZZ
]*x
[n
][ZZ
];
554 /* compute final boxlengths */
555 for (d
=0; d
<DIM
; d
++) {
556 box
[d
][XX
] = mu
[XX
][XX
]*box
[d
][XX
]+mu
[YY
][XX
]*box
[d
][YY
]+mu
[ZZ
][XX
]*box
[d
][ZZ
];
557 box
[d
][YY
] = mu
[YY
][YY
]*box
[d
][YY
]+mu
[ZZ
][YY
]*box
[d
][ZZ
];
558 box
[d
][ZZ
] = mu
[ZZ
][ZZ
]*box
[d
][ZZ
];
561 preserve_box_shape(ir
,box_rel
,box
);
563 /* (un)shifting should NOT be done after this,
564 * since the box vectors might have changed
566 inc_nrnb(nrnb
,eNR_PCOUPL
,nr_atoms
);
569 void berendsen_tcoupl(t_grpopts
*opts
,gmx_ekindata_t
*ekind
,real dt
)
575 for(i
=0; (i
<opts
->ngtc
); i
++) {
576 T
= ekind
->tcstat
[i
].Th
;
578 if ((opts
->tau_t
[i
] > 0) && (T
> 0.0)) {
580 reft
= max(0.0,opts
->ref_t
[i
]);
581 lll
= sqrt(1.0 + (dt
/opts
->tau_t
[i
])*(reft
/T
-1.0));
582 ekind
->tcstat
[i
].lambda
= max(min(lll
,1.25),0.8);
585 ekind
->tcstat
[i
].lambda
= 1.0;
589 fprintf(debug
,"TC: group %d: T: %g, Lambda: %g\n",
590 i
,T
,ekind
->tcstat
[i
].lambda
);
594 void nosehoover_tcoupl(t_grpopts
*opts
,gmx_ekindata_t
*ekind
,real dt
,
595 double xi
[],double vxi
[], t_extmass
*MassQ
)
600 /* note that this routine does not include Nose-hoover chains yet. Should be easy to add. */
602 for(i
=0; (i
<opts
->ngtc
); i
++) {
603 reft
= max(0.0,opts
->ref_t
[i
]);
605 vxi
[i
] += dt
*MassQ
->Qinv
[i
]*(ekind
->tcstat
[i
].Th
- reft
);
606 xi
[i
] += dt
*(oldvxi
+ vxi
[i
])*0.5;
610 t_state
*init_bufstate(const t_state
*template_state
)
613 int nc
= template_state
->nhchainlength
;
615 snew(state
->nosehoover_xi
,nc
*template_state
->ngtc
);
616 snew(state
->nosehoover_vxi
,nc
*template_state
->ngtc
);
617 snew(state
->therm_integral
,template_state
->ngtc
);
618 snew(state
->nhpres_xi
,nc
*template_state
->nnhpres
);
619 snew(state
->nhpres_vxi
,nc
*template_state
->nnhpres
);
624 void destroy_bufstate(t_state
*state
)
628 sfree(state
->nosehoover_xi
);
629 sfree(state
->nosehoover_vxi
);
630 sfree(state
->therm_integral
);
631 sfree(state
->nhpres_xi
);
632 sfree(state
->nhpres_vxi
);
636 void trotter_update(t_inputrec
*ir
,gmx_large_int_t step
, gmx_ekindata_t
*ekind
,
637 gmx_enerdata_t
*enerd
, t_state
*state
,
638 tensor vir
, t_mdatoms
*md
,
639 t_extmass
*MassQ
, int *trotter_seq
)
642 int n
,i
,j
,d
,ntgrp
,ngtc
,gc
=0;
643 t_grp_tcstat
*tcstat
;
645 real ecorr
,pcorr
,dvdlcorr
;
646 real bmass
,qmass
,reft
,kT
,dt
,nd
;
647 tensor dumpres
,dumvir
;
648 double *scalefac
,dtc
;
652 bCouple
= (ir
->nstcalcenergy
== 1 ||
653 do_per_step(step
+ir
->nstcalcenergy
,ir
->nstcalcenergy
));
655 /* signal we are returning if nothing is going to be done in this routine */
656 if ((trotter_seq
[0] == etrtSKIPALL
) || !(bCouple
))
661 dtc
= ir
->nstcalcenergy
*ir
->delta_t
;
662 opts
= &(ir
->opts
); /* just for ease of referencing */
664 snew(scalefac
,opts
->ngtc
);
669 /* execute the series of trotter updates specified in the trotterpart array */
671 for (i
=0;i
<NTROTTERPARTS
;i
++){
672 /* allow for doubled intgrators by doubling dt instead of making 2 calls */
673 if ((trotter_seq
[i
] == etrtBAROV2
) || (trotter_seq
[i
] == etrtBARONHC2
) || (trotter_seq
[i
] == etrtNHC2
))
682 switch (trotter_seq
[i
])
686 boxv_trotter(ir
,&(state
->veta
),dt
,state
->box
,ekind
,vir
,
687 enerd
->term
[F_PDISPCORR
],enerd
->term
[F_DISPCORR
],MassQ
);
691 NHC_trotter(opts
,state
->nnhpres
,ekind
,dt
,state
->nhpres_xi
,
692 state
->nhpres_vxi
,NULL
,&(state
->veta
),MassQ
,FALSE
);
696 NHC_trotter(opts
,opts
->ngtc
,ekind
,dt
,state
->nosehoover_xi
,
697 state
->nosehoover_vxi
,scalefac
,NULL
,MassQ
,(ir
->eI
==eiVV
));
698 /* need to rescale the kinetic energies and velocities here. Could
699 scale the velocities later, but we need them scaled in order to
700 produce the correct outputs, so we'll scale them here. */
702 for (i
=0; i
<ngtc
;i
++)
704 tcstat
= &ekind
->tcstat
[i
];
705 tcstat
->vscale_nhc
= scalefac
[i
];
706 tcstat
->ekinscaleh_nhc
*= (scalefac
[i
]*scalefac
[i
]);
707 tcstat
->ekinscalef_nhc
*= (scalefac
[i
]*scalefac
[i
]);
709 /* now that we've scaled the groupwise velocities, we can add them up to get the total */
710 /* but do we actually need the total? */
712 /* modify the velocities as well */
713 for (n
=md
->start
;n
<md
->start
+md
->homenr
;n
++)
721 state
->v
[n
][d
] *= scalefac
[gc
];
728 sumv
[d
] += (state
->v
[n
][d
])/md
->invmass
[n
];
737 /* check for conserved momentum -- worth looking at this again eventually, but not working right now.*/
745 consk
[d
] = sumv
[d
]*exp((1 + 1.0/opts
->nrdf
[0])*((1.0/DIM
)*log(det(state
->box
)/state
->vol0
)) + state
->nosehoover_xi
[0]);
747 fprintf(debug
,"Conserved kappa: %15.8f %15.8f %15.8f\n",consk
[0],consk
[1],consk
[2]);
754 int **init_npt_vars(t_inputrec
*ir
, t_state
*state
, t_extmass
*MassQ
, bool bTrotter
)
756 int n
,i
,j
,d
,ntgrp
,ngtc
,nnhpres
,nh
,gc
=0;
757 t_grp_tcstat
*tcstat
;
759 real ecorr
,pcorr
,dvdlcorr
;
760 real bmass
,qmass
,reft
,kT
,dt
,ndj
,nd
;
761 tensor dumpres
,dumvir
;
764 opts
= &(ir
->opts
); /* just for ease of referencing */
766 nnhpres
= state
->nnhpres
;
767 nh
= state
->nhchainlength
;
769 if (ir
->eI
== eiMD
) {
770 snew(MassQ
->Qinv
,ngtc
);
771 for(i
=0; (i
<ngtc
); i
++)
773 if ((opts
->tau_t
[i
] > 0) && (opts
->ref_t
[i
] > 0))
775 MassQ
->Qinv
[i
]=1.0/(sqr(opts
->tau_t
[i
]/M_2PI
)*opts
->ref_t
[i
]);
783 else if (EI_VV(ir
->eI
))
785 /* Set pressure variables */
787 if (state
->vol0
== 0)
789 state
->vol0
= det(state
->box
); /* because we start by defining a fixed compressibility,
790 we need the volume at this compressibility to solve the problem */
793 /* units are nm^3 * ns^2 / (nm^3 * bar / kJ/mol) = kJ/mol */
794 /* Investigate this more -- is this the right mass to make it? */
795 MassQ
->Winv
= (PRESFAC
*trace(ir
->compress
)*BOLTZ
*opts
->ref_t
[0])/(DIM
*state
->vol0
*sqr(ir
->tau_p
/M_2PI
));
796 /* An alternate mass definition, from Tuckerman et al. */
797 /* MassQ->Winv = 1.0/(sqr(ir->tau_p/M_2PI)*(opts->nrdf[0]+DIM)*BOLTZ*opts->ref_t[0]); */
802 MassQ
->Winvm
[d
][n
]= PRESFAC
*ir
->compress
[d
][n
]/(state
->vol0
*sqr(ir
->tau_p
/M_2PI
));
803 /* not clear this is correct yet for the anisotropic case*/
806 /* Allocate space for thermostat variables */
807 snew(MassQ
->Qinv
,ngtc
*nh
);
809 /* now, set temperature variables */
810 for(i
=0; i
<ngtc
; i
++)
812 if ((opts
->tau_t
[i
] > 0) && (opts
->ref_t
[i
] > 0))
814 reft
= max(0.0,opts
->ref_t
[i
]);
827 MassQ
->Qinv
[i
*nh
+j
] = 1.0/(sqr(opts
->tau_t
[i
]/M_2PI
)*ndj
*kT
);
835 MassQ
->Qinv
[i
*nh
+j
] = 0.0;
841 /* first, initialize clear all the trotter calls */
842 snew(trotter_seq
,NTROTTERCALLS
);
843 for (i
=0;i
<NTROTTERCALLS
;i
++)
845 snew(trotter_seq
[i
],NTROTTERPARTS
);
846 for (j
=0;j
<NTROTTERPARTS
;j
++) {
847 trotter_seq
[i
][j
] = etrtNONE
;
849 trotter_seq
[i
][0] = etrtSKIPALL
;
854 /* no trotter calls, so we never use the values in the array.
855 * We access them (so we need to define them, but ignore
861 /* compute the kinetic energy by using the half step velocities or
862 * the kinetic energies, depending on the order of the trotter calls */
866 if (IR_NPT_TROTTER(ir
))
868 /* This is the complicated version - there are 4 possible calls, depending on ordering.
869 We start with the initial one. */
870 /* first, a round that estimates veta. */
871 trotter_seq
[0][0] = etrtBAROV
;
873 /* trotter_seq[1] is etrtNHC for 1/2 step velocities - leave zero */
875 /* The first half trotter update */
876 trotter_seq
[2][0] = etrtBAROV
;
877 trotter_seq
[2][1] = etrtNHC
;
878 trotter_seq
[2][2] = etrtBARONHC
;
880 /* The second half trotter update */
881 trotter_seq
[3][0] = etrtBARONHC
;
882 trotter_seq
[3][1] = etrtNHC
;
883 trotter_seq
[3][2] = etrtBAROV
;
885 /* trotter_seq[4] is etrtNHC for second 1/2 step velocities - leave zero */
890 if (IR_NVT_TROTTER(ir
))
892 /* This is the easy version - there are only two calls, both the same.
893 Otherwise, even easier -- no calls */
894 trotter_seq
[2][0] = etrtNHC
;
895 trotter_seq
[3][0] = etrtNHC
;
898 } else if (ir
->eI
==eiVVAK
) {
899 if (IR_NPT_TROTTER(ir
))
901 /* This is the complicated version - there are 4 possible calls, depending on ordering.
902 We start with the initial one. */
903 /* first, a round that estimates veta. */
904 trotter_seq
[0][0] = etrtBAROV
;
906 /* The first half trotter update, part 1 -- double update, because it commutes */
907 trotter_seq
[1][0] = etrtNHC
;
909 /* The first half trotter update, part 2 */
910 trotter_seq
[2][0] = etrtBAROV
;
911 trotter_seq
[2][1] = etrtBARONHC
;
913 /* The second half trotter update, part 1 */
914 trotter_seq
[3][0] = etrtBARONHC
;
915 trotter_seq
[3][1] = etrtBAROV
;
917 /* The second half trotter update -- blank for now */
918 trotter_seq
[4][0] = etrtNHC
;
922 if (IR_NVT_TROTTER(ir
))
924 /* This is the easy version - there is only one call, both the same.
925 Otherwise, even easier -- no calls */
926 trotter_seq
[1][0] = etrtNHC
;
927 trotter_seq
[4][0] = etrtNHC
;
936 bmass
= DIM
*DIM
; /* recommended mass parameters for isotropic barostat */
939 snew(MassQ
->QPinv
,nnhpres
*opts
->nhchainlength
);
941 /* barostat temperature */
942 if ((ir
->tau_p
> 0) && (opts
->ref_t
[0] > 0))
944 reft
= max(0.0,opts
->ref_t
[0]);
946 for (i
=0;i
<nnhpres
;i
++) {
956 MassQ
->QPinv
[i
*opts
->nhchainlength
+j
] = 1.0/(sqr(opts
->tau_t
[0]/M_2PI
)*qmass
*kT
);
962 for (i
=0;i
<nnhpres
;i
++) {
965 MassQ
->QPinv
[i
*nh
+j
]=0.0;
972 real
NPT_energy(t_inputrec
*ir
, t_state
*state
, t_extmass
*MassQ
)
974 int i
,j
,nd
,ndj
,bmass
,qmass
,ngtcall
;
975 real ener_npt
,reft
,eta
,kT
,tau
;
979 int nh
= state
->nhchainlength
;
983 /* now we compute the contribution of the pressure to the conserved quantity*/
985 if (ir
->epc
==epcMTTK
)
987 /* find the volume, and the kinetic energy of the volume */
992 /* contribution from the pressure momenenta */
993 ener_npt
+= 0.5*sqr(state
->veta
)/MassQ
->Winv
;
995 /* contribution from the PV term */
996 vol
= det(state
->box
);
997 ener_npt
+= vol
*trace(ir
->ref_p
)/(DIM
*PRESFAC
);
1000 case epctANISOTROPIC
:
1004 case epctSURFACETENSION
:
1007 case epctSEMIISOTROPIC
:
1015 if (IR_NPT_TROTTER(ir
))
1017 /* add the energy from the barostat thermostat chain */
1018 for (i
=0;i
<state
->nnhpres
;i
++) {
1020 /* note -- assumes only one degree of freedom that is thermostatted in barostat */
1021 ivxi
= &state
->nhpres_vxi
[i
*nh
];
1022 ixi
= &state
->nhpres_xi
[i
*nh
];
1023 iQinv
= &(MassQ
->QPinv
[i
*nh
]);
1024 reft
= max(ir
->opts
.ref_t
[0],0); /* using 'System' temperature */
1031 ener_npt
+= 0.5*sqr(ivxi
[j
])/iQinv
[j
];
1032 /* contribution from the thermal variable of the NH chain */
1033 ener_npt
+= ixi
[j
]*kT
;
1037 fprintf(debug
,"P-T-group: %10d Chain %4d ThermV: %15.8f ThermX: %15.8f",i
,j
,ivxi
[j
],ixi
[j
]);
1045 for(i
=0; i
<ir
->opts
.ngtc
; i
++)
1047 ixi
= &state
->nosehoover_xi
[i
*nh
];
1048 ivxi
= &state
->nosehoover_vxi
[i
*nh
];
1049 iQinv
= &(MassQ
->Qinv
[i
*nh
]);
1051 nd
= ir
->opts
.nrdf
[i
];
1052 reft
= max(ir
->opts
.ref_t
[i
],0);
1057 if (IR_NVT_TROTTER(ir
))
1059 /* contribution from the thermal momenta of the NH chain */
1063 ener_npt
+= 0.5*sqr(ivxi
[j
])/iQinv
[j
];
1064 /* contribution from the thermal variable of the NH chain */
1072 ener_npt
+= ndj
*ixi
[j
]*kT
;
1076 else /* Other non Trotter temperature NH control -- no chains yet. */
1078 ener_npt
+= 0.5*BOLTZ
*nd
*sqr(ivxi
[0])/iQinv
[0];
1079 ener_npt
+= nd
*ixi
[0]*kT
;
1087 static real
vrescale_gamdev(int ia
, gmx_rng_t rng
)
1088 /* Gamma distribution, adapted from numerical recipes */
1091 real am
,e
,s
,v1
,v2
,x
,y
;
1095 for(j
=1; j
<=ia
; j
++) {
1096 x
*= gmx_rng_uniform_real(rng
);
1103 v1
= gmx_rng_uniform_real(rng
);
1104 v2
= 2.0*gmx_rng_uniform_real(rng
)-1.0;
1105 } while (v1
*v1
+ v2
*v2
> 1.0);
1108 s
= sqrt(2.0*am
+ 1.0);
1111 e
= (1.0 + y
*y
)*exp(am
*log(x
/am
) - s
*y
);
1112 } while (gmx_rng_uniform_real(rng
) > e
);
1118 static real
vrescale_sumnoises(int nn
,gmx_rng_t rng
)
1121 * Returns the sum of n independent gaussian noises squared
1122 * (i.e. equivalent to summing the square of the return values
1123 * of nn calls to gmx_rng_gaussian_real).xs
1129 } else if (nn
== 1) {
1130 rr
= gmx_rng_gaussian_real(rng
);
1132 } else if (nn
% 2 == 0) {
1133 return 2.0*vrescale_gamdev(nn
/2,rng
);
1135 rr
= gmx_rng_gaussian_real(rng
);
1136 return 2.0*vrescale_gamdev((nn
-1)/2,rng
) + rr
*rr
;
1140 static real
vrescale_resamplekin(real kk
,real sigma
, int ndeg
, real taut
,
1144 * Generates a new value for the kinetic energy,
1145 * according to Bussi et al JCP (2007), Eq. (A7)
1146 * kk: present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
1147 * sigma: target average value of the kinetic energy (ndeg k_b T/2) (in the same units as kk)
1148 * ndeg: number of degrees of freedom of the atoms to be thermalized
1149 * taut: relaxation time of the thermostat, in units of 'how often this routine is called'
1154 factor
= exp(-1.0/taut
);
1158 rr
= gmx_rng_gaussian_real(rng
);
1161 (1.0 - factor
)*(sigma
*(vrescale_sumnoises(ndeg
-1,rng
) + rr
*rr
)/ndeg
- kk
) +
1162 2.0*rr
*sqrt(kk
*sigma
/ndeg
*(1.0 - factor
)*factor
);
1165 void vrescale_tcoupl(t_grpopts
*opts
,gmx_ekindata_t
*ekind
,real dt
,
1166 double therm_integral
[],gmx_rng_t rng
)
1169 real Ek
,Ek_ref1
,Ek_ref
,Ek_new
;
1171 for(i
=0; (i
<opts
->ngtc
); i
++) {
1172 Ek
= trace(ekind
->tcstat
[i
].ekinh
);
1174 if (opts
->tau_t
[i
] >= 0 && opts
->nrdf
[i
] > 0 && Ek
> 0) {
1175 Ek_ref1
= 0.5*opts
->ref_t
[i
]*BOLTZ
;
1176 Ek_ref
= Ek_ref1
*opts
->nrdf
[i
];
1179 vrescale_resamplekin(Ek
,Ek_ref
,opts
->nrdf
[i
],opts
->tau_t
[i
]/dt
,rng
);
1181 /* Analytically Ek_new>=0, but we check for rounding errors */
1183 ekind
->tcstat
[i
].lambda
= 0.0;
1185 ekind
->tcstat
[i
].lambda
= sqrt(Ek_new
/Ek
);
1187 therm_integral
[i
] -= Ek_new
- Ek
;
1190 fprintf(debug
,"TC: group %d: Ekr %g, Ek %g, Ek_new %g, Lambda: %g\n",
1191 i
,Ek_ref
,Ek
,Ek_new
,ekind
->tcstat
[i
].lambda
);
1194 ekind
->tcstat
[i
].lambda
= 1.0;
1199 real
vrescale_energy(t_grpopts
*opts
,double therm_integral
[])
1205 for(i
=0; i
<opts
->ngtc
; i
++) {
1206 ener
+= therm_integral
[i
];
1212 /* set target temperatures if we are annealing */
1213 void update_annealing_target_temp(t_grpopts
*opts
,real t
)
1216 real pert
,thist
=0,x
;
1218 for(i
=0;i
<opts
->ngtc
;i
++) {
1219 npoints
= opts
->anneal_npoints
[i
];
1220 switch (opts
->annealing
[i
]) {
1224 /* calculate time modulo the period */
1225 pert
= opts
->anneal_time
[i
][npoints
-1];
1227 thist
= t
- n
*pert
; /* modulo time */
1228 /* Make sure rounding didn't get us outside the interval */
1229 if (fabs(thist
-pert
) < GMX_REAL_EPS
*100)
1236 gmx_fatal(FARGS
,"Death horror in update_annealing_target_temp (i=%d/%d npoints=%d)",i
,opts
->ngtc
,npoints
);
1238 /* We are doing annealing for this group if we got here,
1239 * and we have the (relative) time as thist.
1240 * calculate target temp */
1242 while ((j
< npoints
-1) && (thist
>(opts
->anneal_time
[i
][j
+1])))
1244 if (j
< npoints
-1) {
1245 /* Found our position between points j and j+1.
1246 * Interpolate: x is the amount from j+1, (1-x) from point j
1247 * First treat possible jumps in temperature as a special case.
1249 if ((opts
->anneal_time
[i
][j
+1]-opts
->anneal_time
[i
][j
]) < GMX_REAL_EPS
*100)
1250 opts
->ref_t
[i
]=opts
->anneal_temp
[i
][j
+1];
1252 x
= ((thist
-opts
->anneal_time
[i
][j
])/
1253 (opts
->anneal_time
[i
][j
+1]-opts
->anneal_time
[i
][j
]));
1254 opts
->ref_t
[i
] = x
*opts
->anneal_temp
[i
][j
+1]+(1-x
)*opts
->anneal_temp
[i
][j
];
1258 opts
->ref_t
[i
] = opts
->anneal_temp
[i
][npoints
-1];