?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Changed the pV term to the correct one, with the reference pressure. Also made it
[gromacs/rigid-bodies.git] / src / mdlib / mdebin.c
blobe2dd371f24b9956f36afdcf94ce694da2f79110c
1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
2 *
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 3.2.0
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
15
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
20 *
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
27 *
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
30 *
31 * For more info, check our website at http://www.gromacs.org
32 *
33 * And Hey:
34 * GROwing Monsters And Cloning Shrimps
35 */
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 #include <string.h>
41 #include "typedefs.h"
42 #include "string2.h"
43 #include "mdebin.h"
44 #include "smalloc.h"
45 #include "physics.h"
46 #include "enxio.h"
47 #include "vec.h"
48 #include "disre.h"
49 #include "main.h"
50 #include "network.h"
51 #include "names.h"
52 #include "orires.h"
53 #include "constr.h"
54 #include "mtop_util.h"
55 #include "xvgr.h"
56 #include "gmxfio.h"
57
58 static const char *conrmsd_nm[] = { "Constr. rmsd", "Constr.2 rmsd" };
59
60 static const char *boxs_nm[] = { "Box-X", "Box-Y", "Box-Z" };
61
62 static const char *tricl_boxs_nm[] = { "Box-XX", "Box-YX", "Box-YY",
63 "Box-ZX", "Box-ZY", "Box-ZZ" };
64
65 static const char *vol_nm[] = { "Volume" };
66
67 static const char *dens_nm[] = {"Density" };
68
69 static const char *pv_nm[] = {"pV" };
70
71 static const char *boxvel_nm[] = {
72 "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
73 "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
74 };
75
76 #define NBOXS asize(boxs_nm)
77 #define NTRICLBOXS asize(tricl_boxs_nm)
78
79 t_mdebin *init_mdebin(ener_file_t fp_ene,
80 const gmx_mtop_t *mtop,
81 const t_inputrec *ir)
82 {
83 const char *ener_nm[F_NRE];
84 static const char *vir_nm[] = {
85 "Vir-XX", "Vir-XY", "Vir-XZ",
86 "Vir-YX", "Vir-YY", "Vir-YZ",
87 "Vir-ZX", "Vir-ZY", "Vir-ZZ"
88 };
89 static const char *sv_nm[] = {
90 "ShakeVir-XX", "ShakeVir-XY", "ShakeVir-XZ",
91 "ShakeVir-YX", "ShakeVir-YY", "ShakeVir-YZ",
92 "ShakeVir-ZX", "ShakeVir-ZY", "ShakeVir-ZZ"
93 };
94 static const char *fv_nm[] = {
95 "ForceVir-XX", "ForceVir-XY", "ForceVir-XZ",
96 "ForceVir-YX", "ForceVir-YY", "ForceVir-YZ",
97 "ForceVir-ZX", "ForceVir-ZY", "ForceVir-ZZ"
98 };
99 static const char *pres_nm[] = {
100 "Pres-XX","Pres-XY","Pres-XZ",
101 "Pres-YX","Pres-YY","Pres-YZ",
102 "Pres-ZX","Pres-ZY","Pres-ZZ"
103 };
104 static const char *surft_nm[] = {
105 "#Surf*SurfTen"
106 };
107 static const char *mu_nm[] = {
108 "Mu-X", "Mu-Y", "Mu-Z"
109 };
110 static const char *vcos_nm[] = {
111 "2CosZ*Vel-X"
112 };
113 static const char *visc_nm[] = {
114 "1/Viscosity"
115 };
116 char **grpnms;
117 const gmx_groups_t *groups;
118 char **gnm;
119 char buf[256];
120 t_mdebin *md;
121 int i,j,ni,nj,n,k,kk,ncon,nset;
122 bool bBHAM,b14;
123
124 snew(md,1);
125
126 groups = &mtop->groups;
127
128 bBHAM = (mtop->ffparams.functype[0] == F_BHAM);
129 b14 = (gmx_mtop_ftype_count(mtop,F_LJ14) > 0 ||
130 gmx_mtop_ftype_count(mtop,F_LJC14_Q) > 0);
131
132 ncon = gmx_mtop_ftype_count(mtop,F_CONSTR);
133 nset = gmx_mtop_ftype_count(mtop,F_SETTLE);
134 md->bConstr = (ncon > 0 || nset > 0);
135 md->bConstrVir = FALSE;
136 if (md->bConstr) {
137 if (ncon > 0 && ir->eConstrAlg == econtLINCS) {
138 if (ir->eI == eiSD2)
139 md->nCrmsd = 2;
140 else
141 md->nCrmsd = 1;
142 }
143 md->bConstrVir = (getenv("GMX_CONSTRAINTVIR") != NULL);
144 } else {
145 md->nCrmsd = 0;
146 }
147
148 /* Energy monitoring */
149 for(i=0;i<egNR;i++)
150 {
151 md->bEInd[i]=FALSE;
152 }
153
154
155 for(i=0; i<F_NRE; i++) {
156 md->bEner[i] = FALSE;
157 if (i == F_LJ)
158 md->bEner[i] = !bBHAM;
159 else if (i == F_BHAM)
160 md->bEner[i] = bBHAM;
161 else if (i == F_EQM)
162 md->bEner[i] = ir->bQMMM;
163 else if (i == F_COUL_LR)
164 md->bEner[i] = (ir->rcoulomb > ir->rlist);
165 else if (i == F_LJ_LR)
166 md->bEner[i] = (!bBHAM && ir->rvdw > ir->rlist);
167 else if (i == F_BHAM_LR)
168 md->bEner[i] = (bBHAM && ir->rvdw > ir->rlist);
169 else if (i == F_RF_EXCL)
170 md->bEner[i] = (EEL_RF(ir->coulombtype) && ir->coulombtype != eelRF_NEC);
171 else if (i == F_COUL_RECIP)
172 md->bEner[i] = EEL_FULL(ir->coulombtype);
173 else if (i == F_LJ14)
174 md->bEner[i] = b14;
175 else if (i == F_COUL14)
176 md->bEner[i] = b14;
177 else if (i == F_LJC14_Q || i == F_LJC_PAIRS_NB)
178 md->bEner[i] = FALSE;
179 else if ((i == F_DVDL) || (i == F_DKDL))
180 md->bEner[i] = (ir->efep != efepNO);
181 else if (i == F_DHDL_CON)
182 md->bEner[i] = (ir->efep != efepNO && md->bConstr);
183 else if ((interaction_function[i].flags & IF_VSITE) ||
184 (i == F_CONSTR) || (i == F_CONSTRNC) || (i == F_SETTLE))
185 md->bEner[i] = FALSE;
186 else if ((i == F_COUL_SR) || (i == F_EPOT) || (i == F_PRES) || (i==F_EQM))
187 md->bEner[i] = TRUE;
188 else if ((i == F_ETOT) || (i == F_EKIN) || (i == F_TEMP))
189 md->bEner[i] = EI_DYNAMICS(ir->eI);
190 else if (i == F_DISPCORR || i == F_PDISPCORR)
191 md->bEner[i] = (ir->eDispCorr != edispcNO);
192 else if (i == F_DISRESVIOL)
193 md->bEner[i] = (gmx_mtop_ftype_count(mtop,F_DISRES) > 0);
194 else if (i == F_ORIRESDEV)
195 md->bEner[i] = (gmx_mtop_ftype_count(mtop,F_ORIRES) > 0);
196 else if (i == F_CONNBONDS)
197 md->bEner[i] = FALSE;
198 else if (i == F_COM_PULL)
199 md->bEner[i] = (ir->ePull == epullUMBRELLA || ir->ePull == epullCONST_F);
200 else if (i == F_ECONSERVED)
201 md->bEner[i] = ((ir->etc == etcNOSEHOOVER || ir->etc == etcVRESCALE) &&
202 ir->epc == epcNO);
203 else
204 md->bEner[i] = (gmx_mtop_ftype_count(mtop,i) > 0);
205 }
206
207 md->f_nre=0;
208 for(i=0; i<F_NRE; i++)
209 {
210 if (md->bEner[i])
211 {
212 /* FIXME: The constness should not be cast away */
213 /*ener_nm[f_nre]=(char *)interaction_function[i].longname;*/
214 ener_nm[md->f_nre]=interaction_function[i].longname;
215 md->f_nre++;
216 }
217 }
218
219 md->epc = ir->epc;
220 md->bTricl = TRICLINIC(ir->compress) || TRICLINIC(ir->deform);
221 md->bDynBox = DYNAMIC_BOX(*ir);
222 md->etc = ir->etc;
223
224 md->ebin = mk_ebin();
225 /* Pass NULL for unit to let get_ebin_space determine the units
226 * for interaction_function[i].longname
227 */
228 md->ie = get_ebin_space(md->ebin,md->f_nre,ener_nm,NULL);
229 if (md->nCrmsd)
230 {
231 /* This should be called directly after the call for md->ie,
232 * such that md->iconrmsd follows directly in the list.
233 */
234 md->iconrmsd = get_ebin_space(md->ebin,md->nCrmsd,conrmsd_nm,"");
235 }
236 if (md->bDynBox)
237 {
238 md->ib = get_ebin_space(md->ebin, md->bTricl ? NTRICLBOXS :
239 NBOXS, md->bTricl ? tricl_boxs_nm : boxs_nm,
240 unit_length);
241 md->ivol = get_ebin_space(md->ebin, 1, vol_nm, unit_volume);
242 md->idens = get_ebin_space(md->ebin, 1, dens_nm, unit_density_SI);
243 md->ipv = get_ebin_space(md->ebin, 1, pv_nm, unit_energy);
244 }
245 if (md->bConstrVir)
246 {
247 md->isvir = get_ebin_space(md->ebin,asize(sv_nm),sv_nm,unit_energy);
248 md->ifvir = get_ebin_space(md->ebin,asize(fv_nm),fv_nm,unit_energy);
249 }
250 md->ivir = get_ebin_space(md->ebin,asize(vir_nm),vir_nm,unit_energy);
251 md->ipres = get_ebin_space(md->ebin,asize(pres_nm),pres_nm,unit_pres_bar);
252 md->isurft = get_ebin_space(md->ebin,asize(surft_nm),surft_nm,
253 unit_surft_bar);
254 if (md->epc == epcPARRINELLORAHMAN)
255 {
256 md->ipc = get_ebin_space(md->ebin,md->bTricl ? 6 : 3,
257 boxvel_nm,unit_vel);
258 }
259 md->imu = get_ebin_space(md->ebin,asize(mu_nm),mu_nm,unit_dipole_D);
260 if (ir->cos_accel != 0)
261 {
262 md->ivcos = get_ebin_space(md->ebin,asize(vcos_nm),vcos_nm,unit_vel);
263 md->ivisc = get_ebin_space(md->ebin,asize(visc_nm),visc_nm,
264 unit_invvisc_SI);
265 }
266 if (ir->rcoulomb > ir->rlist)
267 {
268 md->bEInd[egCOULLR] = TRUE;
269 }
270 if (!bBHAM)
271 {
272 if (ir->rvdw > ir->rlist)
273 {
274 md->bEInd[egLJLR] = TRUE;
275 }
276 }
277 else
278 {
279 md->bEInd[egLJSR] = FALSE;
280 md->bEInd[egBHAMSR] = TRUE;
281 if (ir->rvdw > ir->rlist)
282 {
283 md->bEInd[egBHAMLR] = TRUE;
284 }
285 }
286 if (b14)
287 {
288 md->bEInd[egLJ14] = TRUE;
289 md->bEInd[egCOUL14] = TRUE;
290 }
291 md->nEc=0;
292 for(i=0; (i<egNR); i++)
293 {
294 if (md->bEInd[i])
295 {
296 md->nEc++;
297 }
298 }
299
300 n=groups->grps[egcENER].nr;
301 md->nEg=n;
302 md->nE=(n*(n+1))/2;
303 snew(md->igrp,md->nE);
304 if (md->nE > 1)
305 {
306 n=0;
307 snew(gnm,md->nEc);
308 for(k=0; (k<md->nEc); k++)
309 {
310 snew(gnm[k],STRLEN);
311 }
312 for(i=0; (i<groups->grps[egcENER].nr); i++)
313 {
314 ni=groups->grps[egcENER].nm_ind[i];
315 for(j=i; (j<groups->grps[egcENER].nr); j++)
316 {
317 nj=groups->grps[egcENER].nm_ind[j];
318 for(k=kk=0; (k<egNR); k++)
319 {
320 if (md->bEInd[k])
321 {
322 sprintf(gnm[kk],"%s:%s-%s",egrp_nm[k],
323 *(groups->grpname[ni]),*(groups->grpname[nj]));
324 kk++;
325 }
326 }
327 md->igrp[n]=get_ebin_space(md->ebin,md->nEc,
328 (const char **)gnm,unit_energy);
329 n++;
330 }
331 }
332 for(k=0; (k<md->nEc); k++)
333 {
334 sfree(gnm[k]);
335 }
336 sfree(gnm);
337
338 if (n != md->nE)
339 {
340 gmx_incons("Number of energy terms wrong");
341 }
342 }
343
344 md->nTC=groups->grps[egcTC].nr;
345 snew(md->grpnms,md->nTC);
346 snew(md->tmp_r,md->nTC);
347 snew(md->tmp_v,md->nTC);
348 grpnms = md->grpnms;
349 for(i=0; (i<md->nTC); i++)
350 {
351 ni=groups->grps[egcTC].nm_ind[i];
352 sprintf(buf,"T-%s",*(groups->grpname[ni]));
353 grpnms[i]=strdup(buf);
354 }
355 md->itemp=get_ebin_space(md->ebin,md->nTC,(const char **)grpnms,
356 unit_temp_K);
357 if (md->etc == etcNOSEHOOVER)
358 {
359 for(i=0; (i<md->nTC); i++)
360 {
361 ni=groups->grps[egcTC].nm_ind[i];
362 sprintf(buf,"Xi-%s",*(groups->grpname[ni]));
363 grpnms[i]=strdup(buf);
364 }
365 md->itc=get_ebin_space(md->ebin,md->nTC,(const char **)grpnms,
366 unit_invtime);
367 }
368 else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
369 md->etc == etcVRESCALE)
370 {
371 for(i=0; (i<md->nTC); i++)
372 {
373 ni=groups->grps[egcTC].nm_ind[i];
374 sprintf(buf,"Lamb-%s",*(groups->grpname[ni]));
375 grpnms[i]=strdup(buf);
376 }
377 md->itc=get_ebin_space(md->ebin,md->nTC,(const char **)grpnms,"");
378 }
379 sfree(grpnms);
380
381 md->nU=groups->grps[egcACC].nr;
382 if (md->nU > 1)
383 {
384 snew(grpnms,3*md->nU);
385 for(i=0; (i<md->nU); i++)
386 {
387 ni=groups->grps[egcACC].nm_ind[i];
388 sprintf(buf,"Ux-%s",*(groups->grpname[ni]));
389 grpnms[3*i+XX]=strdup(buf);
390 sprintf(buf,"Uy-%s",*(groups->grpname[ni]));
391 grpnms[3*i+YY]=strdup(buf);
392 sprintf(buf,"Uz-%s",*(groups->grpname[ni]));
393 grpnms[3*i+ZZ]=strdup(buf);
394 }
395 md->iu=get_ebin_space(md->ebin,3*md->nU,(const char **)grpnms,unit_vel);
396 sfree(grpnms);
397 }
398
399 if ( fp_ene )
400 {
401 do_enxnms(fp_ene,&md->ebin->nener,&md->ebin->enm);
402 }
403
404 md->print_grpnms=NULL;
405
406 return md;
407 }
408
409 FILE *open_dhdl(const char *filename,t_inputrec *ir,const output_env_t oenv)
410 {
411 FILE *fp;
412 const char *dhdl="dH/d\\8l\\4",*deltag="\\8D\\4H",*lambda="\\8l\\4";
413 char title[STRLEN],label_x[STRLEN],label_y[STRLEN];
414 int nsets,s;
415 char **setname,buf[STRLEN];
416
417 sprintf(label_x,"%s (%s)","Time",unit_time);
418 if (ir->n_flambda == 0)
419 {
420 sprintf(title,"%s",dhdl);
421 sprintf(label_y,"%s (%s %s)",
422 dhdl,unit_energy,"[\\8l\\4]\\S-1\\N");
423 }
424 else
425 {
426 sprintf(title,"%s, %s",dhdl,deltag);
427 sprintf(label_y,"(%s)",unit_energy);
428 }
429 fp = xvgropen(filename,title,label_x,label_y,oenv);
430
431 if (ir->n_flambda > 0)
432 {
433 /* g_bar has to determine the lambda values used in this simulation
434 * from this xvg legend.
435 */
436 nsets = 1 + ir->n_flambda;
437 snew(setname,nsets);
438 sprintf(buf,"%s %s %g",dhdl,lambda,ir->init_lambda);
439 setname[0] = strdup(buf);
440 for(s=1; s<nsets; s++)
441 {
442 sprintf(buf,"%s %s %g",deltag,lambda,ir->flambda[s-1]);
443 setname[s] = strdup(buf);
444 }
445 xvgr_legend(fp,nsets,setname,oenv);
446
447 for(s=0; s<nsets; s++)
448 {
449 sfree(setname[s]);
450 }
451 sfree(setname);
452 }
453
454 return fp;
455 }
456
457 static void copy_energy(t_mdebin *md, real e[],real ecpy[])
458 {
459 int i,j;
460
461 for(i=j=0; (i<F_NRE); i++)
462 if (md->bEner[i])
463 ecpy[j++] = e[i];
464 if (j != md->f_nre)
465 gmx_incons("Number of energy terms wrong");
466 }
467
468 void upd_mdebin(t_mdebin *md,FILE *fp_dhdl,
469 bool bSum,
470 double time,
471 real tmass,
472 gmx_enerdata_t *enerd,
473 t_state *state,
474 matrix box,
475 tensor svir,
476 tensor fvir,
477 tensor vir,
478 tensor pres,
479 gmx_ekindata_t *ekind,
480 rvec mu_tot,
481 gmx_constr_t constr)
482 {
483 int i,j,k,kk,m,n,gid;
484 real crmsd[2],tmp6[6];
485 real bs[NTRICLBOXS],vol,dens,pv;
486 real eee[egNR];
487 real ecopy[F_NRE];
488 real tmp;
489
490 /* Do NOT use the box in the state variable, but the separate box provided
491 * as an argument. This is because we sometimes need to write the box from
492 * the last timestep to match the trajectory frames.
493 */
494 copy_energy(md, enerd->term,ecopy);
495 add_ebin(md->ebin,md->ie,md->f_nre,ecopy,bSum);
496 if (md->nCrmsd)
497 {
498 crmsd[0] = constr_rmsd(constr,FALSE);
499 if (md->nCrmsd > 1)
500 {
501 crmsd[1] = constr_rmsd(constr,TRUE);
502 }
503 add_ebin(md->ebin,md->iconrmsd,md->nCrmsd,crmsd,FALSE);
504 }
505 if (md->bDynBox)
506 {
507 if(md->bTricl)
508 {
509 bs[0] = box[XX][XX];
510 bs[1] = box[YY][XX];
511 bs[2] = box[YY][YY];
512 bs[3] = box[ZZ][XX];
513 bs[4] = box[ZZ][YY];
514 bs[5] = box[ZZ][ZZ];
515 }
516 else
517 {
518 bs[0] = box[XX][XX];
519 bs[1] = box[YY][YY];
520 bs[2] = box[ZZ][ZZ];
521 }
522 vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
523 dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
524
525 /* This is pV (in kJ/mol). The pressure is the reference pressure */
526 pv = 0;
527 for (i=0;i<DIM;i++) {
528 for (j=0;j<DIM;j++) {
529 pv += box[i][j]*md->ref_p[i][j];
530 }
531 }
532 add_ebin(md->ebin,md->ib ,NBOXS,bs ,bSum);
533 add_ebin(md->ebin,md->ivol ,1 ,&vol ,bSum);
534 add_ebin(md->ebin,md->idens,1 ,&dens,bSum);
535 add_ebin(md->ebin,md->ipv ,1 ,&pv ,bSum);
536 }
537 if (md->bConstrVir)
538 {
539 add_ebin(md->ebin,md->isvir,9,svir[0],bSum);
540 add_ebin(md->ebin,md->ifvir,9,fvir[0],bSum);
541 }
542 add_ebin(md->ebin,md->ivir,9,vir[0],bSum);
543 add_ebin(md->ebin,md->ipres,9,pres[0],bSum);
544 tmp = (pres[ZZ][ZZ]-(pres[XX][XX]+pres[YY][YY])*0.5)*box[ZZ][ZZ];
545 add_ebin(md->ebin,md->isurft,1,&tmp,bSum);
546 if (md->epc == epcPARRINELLORAHMAN)
547 {
548 tmp6[0] = state->boxv[XX][XX];
549 tmp6[1] = state->boxv[YY][YY];
550 tmp6[2] = state->boxv[ZZ][ZZ];
551 tmp6[3] = state->boxv[YY][XX];
552 tmp6[4] = state->boxv[ZZ][XX];
553 tmp6[5] = state->boxv[ZZ][YY];
554 add_ebin(md->ebin,md->ipc,md->bTricl ? 6 : 3,tmp6,bSum);
555 }
556 add_ebin(md->ebin,md->imu,3,mu_tot,bSum);
557 if (ekind && ekind->cosacc.cos_accel != 0)
558 {
559 vol = box[XX][XX]*box[YY][YY]*box[ZZ][ZZ];
560 dens = (tmass*AMU)/(vol*NANO*NANO*NANO);
561 add_ebin(md->ebin,md->ivcos,1,&(ekind->cosacc.vcos),bSum);
562 /* 1/viscosity, unit 1/(kg m^-1 s^-1) */
563 tmp = 1/(ekind->cosacc.cos_accel/(ekind->cosacc.vcos*PICO)
564 *vol*sqr(box[ZZ][ZZ]*NANO/(2*M_PI)));
565 add_ebin(md->ebin,md->ivisc,1,&tmp,bSum);
566 }
567 if (md->nE > 1)
568 {
569 n=0;
570 for(i=0; (i<md->nEg); i++)
571 {
572 for(j=i; (j<md->nEg); j++)
573 {
574 gid=GID(i,j,md->nEg);
575 for(k=kk=0; (k<egNR); k++)
576 {
577 if (md->bEInd[k])
578 {
579 eee[kk++] = enerd->grpp.ener[k][gid];
580 }
581 }
582 add_ebin(md->ebin,md->igrp[n],md->nEc,eee,bSum);
583 n++;
584 }
585 }
586 }
587
588 if (ekind)
589 {
590 for(i=0; (i<md->nTC); i++)
591 {
592 md->tmp_r[i] = ekind->tcstat[i].T;
593 }
594 add_ebin(md->ebin,md->itemp,md->nTC,md->tmp_r,bSum);
595 if (md->etc == etcNOSEHOOVER)
596 {
597 for(i=0; (i<md->nTC); i++)
598 {
599 md->tmp_r[i] = state->nosehoover_xi[i];
600 }
601 add_ebin(md->ebin,md->itc,md->nTC,md->tmp_r,bSum);
602 }
603 else if (md->etc == etcBERENDSEN || md->etc == etcYES ||
604 md->etc == etcVRESCALE)
605 {
606 for(i=0; (i<md->nTC); i++)
607 {
608 md->tmp_r[i] = ekind->tcstat[i].lambda;
609 }
610 add_ebin(md->ebin,md->itc,md->nTC,md->tmp_r,bSum);
611 }
612 }
613
614 if (ekind && md->nU > 1)
615 {
616 for(i=0; (i<md->nU); i++)
617 {
618 copy_rvec(ekind->grpstat[i].u,md->tmp_v[i]);
619 }
620 add_ebin(md->ebin,md->iu,3*md->nU,md->tmp_v[0],bSum);
621 }
622
623 ebin_increase_count(md->ebin,bSum);
624
625 if (fp_dhdl)
626 {
627 fprintf(fp_dhdl,"%.4f %g",
628 time,
629 enerd->term[F_DVDL]+enerd->term[F_DKDL]+enerd->term[F_DHDL_CON]);
630 for(i=1; i<enerd->n_lambda; i++)
631 {
632 fprintf(fp_dhdl," %g",
633 enerd->enerpart_lambda[i]-enerd->enerpart_lambda[0]);
634 }
635 fprintf(fp_dhdl,"\n");
636 }
637 }
638
639 void upd_mdebin_step(t_mdebin *md)
640 {
641 ebin_increase_count(md->ebin,FALSE);
642 }
643
644 static void npr(FILE *log,int n,char c)
645 {
646 for(; (n>0); n--) fprintf(log,"%c",c);
647 }
648
649 static void pprint(FILE *log,const char *s,t_mdebin *md)
650 {
651 char CHAR='#';
652 int slen;
653 char buf1[22],buf2[22];
654
655 slen = strlen(s);
656 fprintf(log,"\t<====== ");
657 npr(log,slen,CHAR);
658 fprintf(log," ==>\n");
659 fprintf(log,"\t<==== %s ====>\n",s);
660 fprintf(log,"\t<== ");
661 npr(log,slen,CHAR);
662 fprintf(log," ======>\n\n");
663
664 fprintf(log,"\tStatistics over %s steps using %s frames\n",
665 gmx_step_str(md->ebin->nsteps_sim,buf1),
666 gmx_step_str(md->ebin->nsum_sim,buf2));
667 fprintf(log,"\n");
668 }
669
670 void print_ebin_header(FILE *log,gmx_large_int_t steps,double time,real lamb)
671 {
672 char buf[22];
673
674 fprintf(log," %12s %12s %12s\n"
675 " %12s %12.5f %12.5f\n\n",
676 "Step","Time","Lambda",gmx_step_str(steps,buf),time,lamb);
677 }
678
679 void print_ebin(ener_file_t fp_ene,bool bEne,bool bDR,bool bOR,
680 FILE *log,
681 gmx_large_int_t step,double time,
682 int mode,bool bCompact,
683 t_mdebin *md,t_fcdata *fcd,
684 gmx_groups_t *groups,t_grpopts *opts)
685 {
686 /*static char **grpnms=NULL;*/
687 char buf[246];
688 int i,j,n,ni,nj,ndr,nor;
689 int nr[enxNR];
690 real *block[enxNR];
691 t_enxframe fr;
692
693 switch (mode)
694 {
695 case eprNORMAL:
696 fr.t = time;
697 fr.step = step;
698 fr.nsteps = md->ebin->nsteps;
699 fr.nsum = md->ebin->nsum;
700 fr.nre = (bEne) ? md->ebin->nener : 0;
701 fr.ener = md->ebin->e;
702 fr.ndisre = bDR ? fcd->disres.npair : 0;
703 fr.disre_rm3tav = fcd->disres.rm3tav;
704 fr.disre_rt = fcd->disres.rt;
705 /* Optional additional blocks */
706 for(i=0; i<enxNR; i++)
707 {
708 nr[i] = 0;
709 }
710 if (fcd->orires.nr > 0 && bOR)
711 {
712 diagonalize_orires_tensors(&(fcd->orires));
713 nr[enxOR] = fcd->orires.nr;
714 block[enxOR] = fcd->orires.otav;
715 nr[enxORI] = (fcd->orires.oinsl != fcd->orires.otav) ?
716 fcd->orires.nr : 0;
717 block[enxORI] = fcd->orires.oinsl;
718 nr[enxORT] = fcd->orires.nex*12;
719 block[enxORT] = fcd->orires.eig;
720 }
721 fr.nblock = 0;
722 for(i=0; i<enxNR; i++)
723 {
724 if (nr[i] > 0)
725 {
726 fr.nblock = i + 1;
727 }
728 }
729 fr.nr = nr;
730 fr.block = block;
731 if (fr.nre || fr.ndisre || fr.nr[enxOR] || fr.nr[enxORI])
732 {
733 do_enx(fp_ene,&fr);
734 gmx_fio_check_file_position(enx_file_pointer(fp_ene));
735 if (fr.nre)
736 {
737 /* We have stored the sums, so reset the sum history */
738 reset_ebin_sums(md->ebin);
739 }
740 }
741 break;
742 case eprAVER:
743 if (log)
744 {
745 pprint(log,"A V E R A G E S",md);
746 }
747 break;
748 case eprRMS:
749 if (log)
750 {
751 pprint(log,"R M S - F L U C T U A T I O N S",md);
752 }
753 break;
754 default:
755 gmx_fatal(FARGS,"Invalid print mode (%d)",mode);
756 }
757
758 if (log)
759 {
760 for(i=0;i<opts->ngtc;i++)
761 {
762 if(opts->annealing[i]!=eannNO)
763 {
764 fprintf(log,"Current ref_t for group %s: %8.1f\n",
765 *(groups->grpname[groups->grps[egcTC].nm_ind[i]]),
766 opts->ref_t[i]);
767 }
768 }
769 if (mode==eprNORMAL && fcd->orires.nr>0)
770 {
771 print_orires_log(log,&(fcd->orires));
772 }
773 fprintf(log," Energies (%s)\n",unit_energy);
774 pr_ebin(log,md->ebin,md->ie,md->f_nre+md->nCrmsd,5,mode,TRUE);
775 fprintf(log,"\n");
776
777 if (!bCompact)
778 {
779 if (md->bDynBox)
780 {
781 pr_ebin(log,md->ebin,md->ib, md->bTricl ? NTRICLBOXS : NBOXS,5,
782 mode,TRUE);
783 fprintf(log,"\n");
784 }
785 if (md->bConstrVir)
786 {
787 fprintf(log," Constraint Virial (%s)\n",unit_energy);
788 pr_ebin(log,md->ebin,md->isvir,9,3,mode,FALSE);
789 fprintf(log,"\n");
790 fprintf(log," Force Virial (%s)\n",unit_energy);
791 pr_ebin(log,md->ebin,md->ifvir,9,3,mode,FALSE);
792 fprintf(log,"\n");
793 }
794 fprintf(log," Total Virial (%s)\n",unit_energy);
795 pr_ebin(log,md->ebin,md->ivir,9,3,mode,FALSE);
796 fprintf(log,"\n");
797 fprintf(log," Pressure (%s)\n",unit_pres_bar);
798 pr_ebin(log,md->ebin,md->ipres,9,3,mode,FALSE);
799 fprintf(log,"\n");
800 fprintf(log," Total Dipole (%s)\n",unit_dipole_D);
801 pr_ebin(log,md->ebin,md->imu,3,3,mode,FALSE);
802 fprintf(log,"\n");
803
804 if (md->nE > 1)
805 {
806 if (md->print_grpnms==NULL)
807 {
808 snew(md->print_grpnms,md->nE);
809 n=0;
810 for(i=0; (i<md->nEg); i++)
811 {
812 ni=groups->grps[egcENER].nm_ind[i];
813 for(j=i; (j<md->nEg); j++)
814 {
815 nj=groups->grps[egcENER].nm_ind[j];
816 sprintf(buf,"%s-%s",*(groups->grpname[ni]),
817 *(groups->grpname[nj]));
818 md->print_grpnms[n++]=strdup(buf);
819 }
820 }
821 }
822 sprintf(buf,"Epot (%s)",unit_energy);
823 fprintf(log,"%15s ",buf);
824 for(i=0; (i<egNR); i++)
825 {
826 if (md->bEInd[i])
827 {
828 fprintf(log,"%12s ",egrp_nm[i]);
829 }
830 }
831 fprintf(log,"\n");
832 for(i=0; (i<md->nE); i++)
833 {
834 fprintf(log,"%15s",md->print_grpnms[i]);
835 pr_ebin(log,md->ebin,md->igrp[i],md->nEc,md->nEc,mode,
836 FALSE);
837 }
838 fprintf(log,"\n");
839 }
840 if (md->nTC > 1)
841 {
842 pr_ebin(log,md->ebin,md->itemp,md->nTC,4,mode,TRUE);
843 fprintf(log,"\n");
844 }
845 if (md->nU > 1)
846 {
847 fprintf(log,"%15s %12s %12s %12s\n",
848 "Group","Ux","Uy","Uz");
849 for(i=0; (i<md->nU); i++)
850 {
851 ni=groups->grps[egcACC].nm_ind[i];
852 fprintf(log,"%15s",*groups->grpname[ni]);
853 pr_ebin(log,md->ebin,md->iu+3*i,3,3,mode,FALSE);
854 }
855 fprintf(log,"\n");
856 }
857 }
858
859 }
860 }
861
862 void
863 init_energyhistory(energyhistory_t * enerhist)
864 {
865 enerhist->nener = 0;
866
867 enerhist->ener_ave = NULL;
868 enerhist->ener_sum = NULL;
869 enerhist->ener_sum_sim = NULL;
870
871 enerhist->nsteps = 0;
872 enerhist->nsum = 0;
873 enerhist->nsteps_sim = 0;
874 enerhist->nsum_sim = 0;
875 }
876
877 void
878 update_energyhistory(energyhistory_t * enerhist,t_mdebin * mdebin)
879 {
880 int i;
881
882 enerhist->nsteps = mdebin->ebin->nsteps;
883 enerhist->nsum = mdebin->ebin->nsum;
884 enerhist->nsteps_sim = mdebin->ebin->nsteps_sim;
885 enerhist->nsum_sim = mdebin->ebin->nsum_sim;
886 enerhist->nener = mdebin->ebin->nener;
887
888 if (mdebin->ebin->nsum > 0)
889 {
890 /* Check if we need to allocate first */
891 if(enerhist->ener_ave == NULL)
892 {
893 snew(enerhist->ener_ave,enerhist->nener);
894 snew(enerhist->ener_sum,enerhist->nener);
895 }
896
897 for(i=0;i<enerhist->nener;i++)
898 {
899 enerhist->ener_ave[i] = mdebin->ebin->e[i].eav;
900 enerhist->ener_sum[i] = mdebin->ebin->e[i].esum;
901 }
902 }
903
904 if (mdebin->ebin->nsum_sim > 0)
905 {
906 /* Check if we need to allocate first */
907 if(enerhist->ener_sum_sim == NULL)
908 {
909 snew(enerhist->ener_sum_sim,enerhist->nener);
910 }
911
912 for(i=0;i<enerhist->nener;i++)
913 {
914 enerhist->ener_sum_sim[i] = mdebin->ebin->e_sim[i].esum;
915 }
916 }
917 }
918
919 void
920 restore_energyhistory_from_state(t_mdebin * mdebin,energyhistory_t * enerhist)
921 {
922 int i;
923
924 if ((enerhist->nsum > 0 || enerhist->nsum_sim > 0) &&
925 mdebin->ebin->nener != enerhist->nener)
926 {
927 gmx_fatal(FARGS,"Mismatch between number of energies in run input (%d) and checkpoint file (%d).",
928 mdebin->ebin->nener,enerhist->nener);
929 }
930
931 mdebin->ebin->nsteps = enerhist->nsteps;
932 mdebin->ebin->nsum = enerhist->nsum;
933 mdebin->ebin->nsteps_sim = enerhist->nsteps_sim;
934 mdebin->ebin->nsum_sim = enerhist->nsum_sim;
935
936 for(i=0; i<mdebin->ebin->nener; i++)
937 {
938 mdebin->ebin->e[i].eav =
939 (enerhist->nsum > 0 ? enerhist->ener_ave[i] : 0);
940 mdebin->ebin->e[i].esum =
941 (enerhist->nsum > 0 ? enerhist->ener_sum[i] : 0);
942 mdebin->ebin->e_sim[i].esum =
943 (enerhist->nsum_sim > 0 ? enerhist->ener_sum_sim[i] : 0);
944 }
945 }