?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Enforced rotation: fixed torque calculation for FLEX potential when using mass-weighting
[gromacs/adressmacs.git] / src / tools / gmx_energy.c
blobcd299645b164f26c58635b1cae90ebae8cfa1193
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 * Green Red Orange Magenta Azure Cyan Skyblue
35 */
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 #include <string.h>
41 #include <math.h>
42 #include <ctype.h>
43
44 #include "typedefs.h"
45 #include "gmx_fatal.h"
46 #include "vec.h"
47 #include "string2.h"
48 #include "smalloc.h"
49 #include "enxio.h"
50 #include "statutil.h"
51 #include "names.h"
52 #include "copyrite.h"
53 #include "macros.h"
54 #include "xvgr.h"
55 #include "gstat.h"
56 #include "physics.h"
57 #include "tpxio.h"
58 #include "viewit.h"
59 #include "mtop_util.h"
60 #include "gmx_ana.h"
61
62
63 static real minthird=-1.0/3.0,minsixth=-1.0/6.0;
64
65 typedef struct {
66 real sum;
67 real sum2;
68 } exactsum_t;
69
70 typedef struct {
71 real *ener;
72 exactsum_t *es;
73 gmx_bool bExactStat;
74 double av;
75 double rmsd;
76 double ee;
77 double slope;
78 } enerdat_t;
79
80 typedef struct {
81 gmx_large_int_t nsteps;
82 gmx_large_int_t npoints;
83 int nframes;
84 int *step;
85 int *steps;
86 int *points;
87 enerdat_t *s;
88 } enerdata_t;
89
90 static double mypow(double x,double y)
91 {
92 if (x > 0)
93 return pow(x,y);
94 else
95 return 0.0;
96 }
97
98 static int *select_it(int nre,char *nm[],int *nset)
99 {
100 gmx_bool *bE;
101 int n,k,j,i;
102 int *set;
103 gmx_bool bVerbose = TRUE;
104
105 if ((getenv("VERBOSE")) != NULL)
106 bVerbose = FALSE;
107
108 fprintf(stderr,"Select the terms you want from the following list\n");
109 fprintf(stderr,"End your selection with 0\n");
110
111 if ( bVerbose ) {
112 for(k=0; (k<nre); ) {
113 for(j=0; (j<4) && (k<nre); j++,k++)
114 fprintf(stderr," %3d=%14s",k+1,nm[k]);
115 fprintf(stderr,"\n");
116 }
117 }
118
119 snew(bE,nre);
120 do {
121 if(1 != scanf("%d",&n))
122 {
123 gmx_fatal(FARGS,"Error reading user input");
124 }
125 if ((n>0) && (n<=nre))
126 bE[n-1]=TRUE;
127 } while (n != 0);
128
129 snew(set,nre);
130 for(i=(*nset)=0; (i<nre); i++)
131 if (bE[i])
132 set[(*nset)++]=i;
133
134 sfree(bE);
135
136 return set;
137 }
138
139 static int strcount(const char *s1,const char *s2)
140 {
141 int n=0;
142 while (s1 && s2 && (toupper(s1[n]) == toupper(s2[n])))
143 n++;
144 return n;
145 }
146
147 static void chomp(char *buf)
148 {
149 int len = strlen(buf);
150
151 while ((len > 0) && (buf[len-1] == '\n')) {
152 buf[len-1] = '\0';
153 len--;
154 }
155 }
156
157 static int *select_by_name(int nre,gmx_enxnm_t *nm,int *nset)
158 {
159 gmx_bool *bE;
160 int n,k,kk,j,i,nmatch,nind,nss;
161 int *set;
162 gmx_bool bEOF,bVerbose = TRUE,bLong=FALSE;
163 char *ptr,buf[STRLEN];
164 const char *fm4="%3d %-14s";
165 const char *fm2="%3d %-34s";
166 char **newnm=NULL;
167
168 if ((getenv("VERBOSE")) != NULL)
169 bVerbose = FALSE;
170
171 fprintf(stderr,"\n");
172 fprintf(stderr,"Select the terms you want from the following list by\n");
173 fprintf(stderr,"selecting either (part of) the name or the number or a combination.\n");
174 fprintf(stderr,"End your selection with an empty line or a zero.\n");
175 fprintf(stderr,"-------------------------------------------------------------------\n");
176
177 snew(newnm,nre);
178 j = 0;
179 for(k=0; k<nre; k++) {
180 newnm[k] = strdup(nm[k].name);
181 /* Insert dashes in all the names */
182 while ((ptr = strchr(newnm[k],' ')) != NULL) {
183 *ptr = '-';
184 }
185 if ( bVerbose ) {
186 if (j == 0) {
187 if (k > 0) {
188 fprintf(stderr,"\n");
189 }
190 bLong = FALSE;
191 for(kk=k; kk<k+4; kk++) {
192 if (kk < nre && strlen(nm[kk].name) > 14) {
193 bLong = TRUE;
194 }
195 }
196 } else {
197 fprintf(stderr," ");
198 }
199 if (!bLong) {
200 fprintf(stderr,fm4,k+1,newnm[k]);
201 j++;
202 if (j == 4) {
203 j = 0;
204 }
205 } else {
206 fprintf(stderr,fm2,k+1,newnm[k]);
207 j++;
208 if (j == 2) {
209 j = 0;
210 }
211 }
212 }
213 }
214 if ( bVerbose ) {
215 fprintf(stderr,"\n\n");
216 }
217
218 snew(bE,nre);
219
220 bEOF = FALSE;
221 while (!bEOF && (fgets2(buf,STRLEN-1,stdin))) {
222 /* Remove newlines */
223 chomp(buf);
224
225 /* Remove spaces */
226 trim(buf);
227
228 /* Empty line means end of input */
229 bEOF = (strlen(buf) == 0);
230 if (!bEOF) {
231 ptr = buf;
232 do {
233 if (!bEOF) {
234 /* First try to read an integer */
235 nss = sscanf(ptr,"%d",&nind);
236 if (nss == 1) {
237 /* Zero means end of input */
238 if (nind == 0) {
239 bEOF = TRUE;
240 } else if ((1<=nind) && (nind<=nre)) {
241 bE[nind-1] = TRUE;
242 } else {
243 fprintf(stderr,"number %d is out of range\n",nind);
244 }
245 }
246 else {
247 /* Now try to read a string */
248 i = strlen(ptr);
249 nmatch = 0;
250 for(nind=0; nind<nre; nind++) {
251 if (gmx_strcasecmp(newnm[nind],ptr) == 0) {
252 bE[nind] = TRUE;
253 nmatch++;
254 }
255 }
256 if (nmatch == 0) {
257 i = strlen(ptr);
258 nmatch = 0;
259 for(nind=0; nind<nre; nind++) {
260 if (gmx_strncasecmp(newnm[nind],ptr,i) == 0) {
261 bE[nind] = TRUE;
262 nmatch++;
263 }
264 }
265 if (nmatch == 0) {
266 fprintf(stderr,"String '%s' does not match anything\n",ptr);
267 }
268 }
269 }
270 }
271 /* Look for the first space, and remove spaces from there */
272 if ((ptr = strchr(ptr,' ')) != NULL) {
273 trim(ptr);
274 }
275 } while (!bEOF && (ptr && (strlen(ptr) > 0)));
276 }
277 }
278
279 snew(set,nre);
280 for(i=(*nset)=0; (i<nre); i++)
281 if (bE[i])
282 set[(*nset)++]=i;
283
284 sfree(bE);
285
286 if (*nset == 0)
287 gmx_fatal(FARGS,"No energy terms selected");
288
289 for(i=0; (i<nre); i++)
290 sfree(newnm[i]);
291 sfree(newnm);
292
293 return set;
294 }
295
296 static void get_orires_parms(const char *topnm,
297 int *nor,int *nex,int **label,real **obs)
298 {
299 gmx_mtop_t mtop;
300 gmx_localtop_t *top;
301 t_inputrec ir;
302 t_iparams *ip;
303 int natoms,i;
304 t_iatom *iatom;
305 int nb;
306 matrix box;
307
308 read_tpx(topnm,&ir,box,&natoms,NULL,NULL,NULL,&mtop);
309 top = gmx_mtop_generate_local_top(&mtop,&ir);
310
311 ip = top->idef.iparams;
312 iatom = top->idef.il[F_ORIRES].iatoms;
313
314 /* Count how many distance restraint there are... */
315 nb = top->idef.il[F_ORIRES].nr;
316 if (nb == 0)
317 gmx_fatal(FARGS,"No orientation restraints in topology!\n");
318
319 *nor = nb/3;
320 *nex = 0;
321 snew(*label,*nor);
322 snew(*obs,*nor);
323 for(i=0; i<nb; i+=3) {
324 (*label)[i/3] = ip[iatom[i]].orires.label;
325 (*obs)[i/3] = ip[iatom[i]].orires.obs;
326 if (ip[iatom[i]].orires.ex >= *nex)
327 *nex = ip[iatom[i]].orires.ex+1;
328 }
329 fprintf(stderr,"Found %d orientation restraints with %d experiments",
330 *nor,*nex);
331 }
332
333 static int get_bounds(const char *topnm,
334 real **bounds,int **index,int **dr_pair,int *npairs,
335 gmx_mtop_t *mtop,gmx_localtop_t **ltop,t_inputrec *ir)
336 {
337 gmx_localtop_t *top;
338 t_functype *functype;
339 t_iparams *ip;
340 int natoms,i,j,k,type,ftype,natom;
341 t_ilist *disres;
342 t_iatom *iatom;
343 real *b;
344 int *ind,*pair;
345 int nb,label1;
346 matrix box;
347
348 read_tpx(topnm,ir,box,&natoms,NULL,NULL,NULL,mtop);
349 snew(*ltop,1);
350 top = gmx_mtop_generate_local_top(mtop,ir);
351 *ltop = top;
352
353 functype = top->idef.functype;
354 ip = top->idef.iparams;
355
356 /* Count how many distance restraint there are... */
357 nb=top->idef.il[F_DISRES].nr;
358 if (nb == 0)
359 gmx_fatal(FARGS,"No distance restraints in topology!\n");
360
361 /* Allocate memory */
362 snew(b,nb);
363 snew(ind,nb);
364 snew(pair,nb+1);
365
366 /* Fill the bound array */
367 nb=0;
368 for(i=0; (i<top->idef.ntypes); i++) {
369 ftype = functype[i];
370 if (ftype == F_DISRES) {
371
372 label1 = ip[i].disres.label;
373 b[nb] = ip[i].disres.up1;
374 ind[nb] = label1;
375 nb++;
376 }
377 }
378 *bounds = b;
379
380 /* Fill the index array */
381 label1 = -1;
382 disres = &(top->idef.il[F_DISRES]);
383 iatom = disres->iatoms;
384 for(i=j=k=0; (i<disres->nr); ) {
385 type = iatom[i];
386 ftype = top->idef.functype[type];
387 natom = interaction_function[ftype].nratoms+1;
388 if (label1 != top->idef.iparams[type].disres.label) {
389 pair[j] = k;
390 label1 = top->idef.iparams[type].disres.label;
391 j ++;
392 }
393 k++;
394 i += natom;
395 }
396 pair[j] = k;
397 *npairs = k;
398 if (j != nb)
399 gmx_incons("get_bounds for distance restraints");
400
401 *index = ind;
402 *dr_pair = pair;
403
404 return nb;
405 }
406
407 static void calc_violations(real rt[],real rav3[],int nb,int index[],
408 real bounds[],real *viol,double *st,double *sa)
409 {
410 const real sixth=1.0/6.0;
411 int i,j;
412 double rsum,rav,sumaver,sumt;
413
414 sumaver = 0;
415 sumt = 0;
416 for(i=0; (i<nb); i++) {
417 rsum = 0.0;
418 rav = 0.0;
419 for(j=index[i]; (j<index[i+1]); j++) {
420 if (viol)
421 viol[j] += mypow(rt[j],-3.0);
422 rav += sqr(rav3[j]);
423 rsum += mypow(rt[j],-6);
424 }
425 rsum = max(0.0,mypow(rsum,-sixth)-bounds[i]);
426 rav = max(0.0,mypow(rav, -sixth)-bounds[i]);
427
428 sumt += rsum;
429 sumaver += rav;
430 }
431 *st = sumt;
432 *sa = sumaver;
433 }
434
435 static void analyse_disre(const char *voutfn, int nframes,
436 real violaver[], real bounds[], int index[],
437 int pair[], int nbounds,
438 const output_env_t oenv)
439 {
440 FILE *vout;
441 double sum,sumt,sumaver;
442 int i,j;
443
444 /* Subtract bounds from distances, to calculate violations */
445 calc_violations(violaver,violaver,
446 nbounds,pair,bounds,NULL,&sumt,&sumaver);
447
448 #ifdef DEBUG
449 fprintf(stdout,"\nSum of violations averaged over simulation: %g nm\n",
450 sumaver);
451 fprintf(stdout,"Largest violation averaged over simulation: %g nm\n\n",
452 sumt);
453 #endif
454 vout=xvgropen(voutfn,"r\\S-3\\N average violations","DR Index","nm",
455 oenv);
456 sum = 0.0;
457 sumt = 0.0;
458 for(i=0; (i<nbounds); i++) {
459 /* Do ensemble averaging */
460 sumaver = 0;
461 for(j=pair[i]; (j<pair[i+1]); j++)
462 sumaver += sqr(violaver[j]/nframes);
463 sumaver = max(0.0,mypow(sumaver,minsixth)-bounds[i]);
464
465 sumt += sumaver;
466 sum = max(sum,sumaver);
467 fprintf(vout,"%10d %10.5e\n",index[i],sumaver);
468 }
469 #ifdef DEBUG
470 for(j=0; (j<dr.ndr); j++)
471 fprintf(vout,"%10d %10.5e\n",j,mypow(violaver[j]/nframes,minthird));
472 #endif
473 ffclose(vout);
474
475 fprintf(stdout,"\nSum of violations averaged over simulation: %g nm\n",
476 sumt);
477 fprintf(stdout,"Largest violation averaged over simulation: %g nm\n\n",sum);
478
479 do_view(oenv,voutfn,"-graphtype bar");
480 }
481
482 static void einstein_visco(const char *fn,const char *fni,int nsets,
483 int nframes,real **sum,
484 real V,real T,int nsteps,double time[],
485 const output_env_t oenv)
486 {
487 FILE *fp0,*fp1;
488 double av[4],avold[4];
489 double fac,dt,di;
490 int i,j,m,nf4;
491
492 if (nframes < 1)
493 return;
494
495 dt = (time[1]-time[0]);
496 nf4 = nframes/4+1;
497
498 for(i=0; i<=nsets; i++)
499 avold[i] = 0;
500 fp0=xvgropen(fni,"Shear viscosity integral",
501 "Time (ps)","(kg m\\S-1\\N s\\S-1\\N ps)",oenv);
502 fp1=xvgropen(fn,"Shear viscosity using Einstein relation",
503 "Time (ps)","(kg m\\S-1\\N s\\S-1\\N)",oenv);
504 for(i=1; i<nf4; i++) {
505 fac = dt*nframes/nsteps;
506 for(m=0; m<=nsets; m++)
507 av[m] = 0;
508 for(j=0; j<nframes-i; j++) {
509 for(m=0; m<nsets; m++) {
510 di = sqr(fac*(sum[m][j+i]-sum[m][j]));
511
512 av[m] += di;
513 av[nsets] += di/nsets;
514 }
515 }
516 /* Convert to SI for the viscosity */
517 fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nframes-i);
518 fprintf(fp0,"%10g",time[i]-time[0]);
519 for(m=0; (m<=nsets); m++) {
520 av[m] = fac*av[m];
521 fprintf(fp0," %10g",av[m]);
522 }
523 fprintf(fp0,"\n");
524 fprintf(fp1,"%10g",0.5*(time[i]+time[i-1])-time[0]);
525 for(m=0; (m<=nsets); m++) {
526 fprintf(fp1," %10g",(av[m]-avold[m])/dt);
527 avold[m] = av[m];
528 }
529 fprintf(fp1,"\n");
530 }
531 ffclose(fp0);
532 ffclose(fp1);
533 }
534
535 typedef struct {
536 gmx_large_int_t np;
537 double sum;
538 double sav;
539 double sav2;
540 } ee_sum_t;
541
542 typedef struct {
543 int b;
544 ee_sum_t sum;
545 gmx_large_int_t nst;
546 gmx_large_int_t nst_min;
547 } ener_ee_t;
548
549 static void clear_ee_sum(ee_sum_t *ees)
550 {
551 ees->sav = 0;
552 ees->sav2 = 0;
553 ees->np = 0;
554 ees->sum = 0;
555 }
556
557 static void add_ee_sum(ee_sum_t *ees,double sum,int np)
558 {
559 ees->np += np;
560 ees->sum += sum;
561 }
562
563 static void add_ee_av(ee_sum_t *ees)
564 {
565 double av;
566
567 av = ees->sum/ees->np;
568 ees->sav += av;
569 ees->sav2 += av*av;
570 ees->np = 0;
571 ees->sum = 0;
572 }
573
574 static double calc_ee2(int nb,ee_sum_t *ees)
575 {
576 return (ees->sav2/nb - dsqr(ees->sav/nb))/(nb - 1);
577 }
578
579 static void set_ee_av(ener_ee_t *eee)
580 {
581 if (debug)
582 {
583 char buf[STEPSTRSIZE];
584 fprintf(debug,"Storing average for err.est.: %s steps\n",
585 gmx_step_str(eee->nst,buf));
586 }
587 add_ee_av(&eee->sum);
588 eee->b++;
589 if (eee->b == 1 || eee->nst < eee->nst_min)
590 {
591 eee->nst_min = eee->nst;
592 }
593 eee->nst = 0;
594 }
595
596 static void calc_averages(int nset,enerdata_t *edat,int nbmin,int nbmax)
597 {
598 int nb,i,f,nee;
599 double sum,sum2,sump,see2;
600 gmx_large_int_t steps,np,p,bound_nb;
601 enerdat_t *ed;
602 exactsum_t *es;
603 gmx_bool bAllZero;
604 double x,sx,sy,sxx,sxy;
605 ener_ee_t *eee;
606
607 /* Check if we have exact statistics over all points */
608 for(i=0; i<nset; i++)
609 {
610 ed = &edat->s[i];
611 ed->bExactStat = FALSE;
612 if (edat->npoints > 0)
613 {
614 /* All energy file sum entries 0 signals no exact sums.
615 * But if all energy values are 0, we still have exact sums.
616 */
617 bAllZero = TRUE;
618 for(f=0; f<edat->nframes && !ed->bExactStat; f++)
619 {
620 if (ed->ener[i] != 0)
621 {
622 bAllZero = FALSE;
623 }
624 ed->bExactStat = (ed->es[f].sum != 0);
625 }
626 if (bAllZero)
627 {
628 ed->bExactStat = TRUE;
629 }
630 }
631 }
632
633 snew(eee,nbmax+1);
634 for(i=0; i<nset; i++)
635 {
636 ed = &edat->s[i];
637
638 sum = 0;
639 sum2 = 0;
640 np = 0;
641 sx = 0;
642 sy = 0;
643 sxx = 0;
644 sxy = 0;
645 for(nb=nbmin; nb<=nbmax; nb++)
646 {
647 eee[nb].b = 0;
648 clear_ee_sum(&eee[nb].sum);
649 eee[nb].nst = 0;
650 eee[nb].nst_min = 0;
651 }
652 for(f=0; f<edat->nframes; f++)
653 {
654 es = &ed->es[f];
655
656 if (ed->bExactStat)
657 {
658 /* Add the sum and the sum of variances to the totals. */
659 p = edat->points[f];
660 sump = es->sum;
661 sum2 += es->sum2;
662 if (np > 0)
663 {
664 sum2 += dsqr(sum/np - (sum + es->sum)/(np + p))
665 *np*(np + p)/p;
666 }
667 }
668 else
669 {
670 /* Add a single value to the sum and sum of squares. */
671 p = 1;
672 sump = ed->ener[f];
673 sum2 += dsqr(sump);
674 }
675
676 /* sum has to be increased after sum2 */
677 np += p;
678 sum += sump;
679
680 /* For the linear regression use variance 1/p.
681 * Note that sump is the sum, not the average, so we don't need p*.
682 */
683 x = edat->step[f] - 0.5*(edat->steps[f] - 1);
684 sx += p*x;
685 sy += sump;
686 sxx += p*x*x;
687 sxy += x*sump;
688
689 for(nb=nbmin; nb<=nbmax; nb++)
690 {
691 /* Check if the current end step is closer to the desired
692 * block boundary than the next end step.
693 */
694 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
695 if (eee[nb].nst > 0 &&
696 bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb)
697 {
698 set_ee_av(&eee[nb]);
699 }
700 if (f == 0)
701 {
702 eee[nb].nst = 1;
703 }
704 else
705 {
706 eee[nb].nst += edat->step[f] - edat->step[f-1];
707 }
708 if (ed->bExactStat)
709 {
710 add_ee_sum(&eee[nb].sum,es->sum,edat->points[f]);
711 }
712 else
713 {
714 add_ee_sum(&eee[nb].sum,edat->s[i].ener[f],1);
715 }
716 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
717 if (edat->step[f]*nb >= bound_nb)
718 {
719 set_ee_av(&eee[nb]);
720 }
721 }
722 }
723
724 edat->s[i].av = sum/np;
725 if (ed->bExactStat)
726 {
727 edat->s[i].rmsd = sqrt(sum2/np);
728 }
729 else
730 {
731 edat->s[i].rmsd = sqrt(sum2/np - dsqr(edat->s[i].av));
732 }
733
734 if (edat->nframes > 1)
735 {
736 edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx);
737 }
738 else
739 {
740 edat->s[i].slope = 0;
741 }
742
743 nee = 0;
744 see2 = 0;
745 for(nb=nbmin; nb<=nbmax; nb++)
746 {
747 /* Check if we actually got nb blocks and if the smallest
748 * block is not shorter than 80% of the average.
749 */
750 if (debug)
751 {
752 char buf1[STEPSTRSIZE],buf2[STEPSTRSIZE];
753 fprintf(debug,"Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
754 nb,eee[nb].b,
755 gmx_step_str(eee[nb].nst_min,buf1),
756 gmx_step_str(edat->nsteps,buf2));
757 }
758 if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps)
759 {
760 see2 += calc_ee2(nb,&eee[nb].sum);
761 nee++;
762 }
763 }
764 if (nee > 0)
765 {
766 edat->s[i].ee = sqrt(see2/nee);
767 }
768 else
769 {
770 edat->s[i].ee = -1;
771 }
772 }
773 sfree(eee);
774 }
775
776 static enerdata_t *calc_sum(int nset,enerdata_t *edat,int nbmin,int nbmax)
777 {
778 enerdata_t *esum;
779 enerdat_t *s;
780 int f,i;
781 double sum;
782
783 snew(esum,1);
784 *esum = *edat;
785 snew(esum->s,1);
786 s = &esum->s[0];
787 snew(s->ener,esum->nframes);
788 snew(s->es ,esum->nframes);
789
790 s->bExactStat = TRUE;
791 s->slope = 0;
792 for(i=0; i<nset; i++)
793 {
794 if (!edat->s[i].bExactStat)
795 {
796 s->bExactStat = FALSE;
797 }
798 s->slope += edat->s[i].slope;
799 }
800
801 for(f=0; f<edat->nframes; f++)
802 {
803 sum = 0;
804 for(i=0; i<nset; i++)
805 {
806 sum += edat->s[i].ener[f];
807 }
808 s->ener[f] = sum;
809 sum = 0;
810 for(i=0; i<nset; i++)
811 {
812 sum += edat->s[i].es[f].sum;
813 }
814 s->es[f].sum = sum;
815 s->es[f].sum2 = 0;
816 }
817
818 calc_averages(1,esum,nbmin,nbmax);
819
820 return esum;
821 }
822
823 static char *ee_pr(double ee,char *buf)
824 {
825 char tmp[100];
826 double rnd;
827
828 if (ee < 0)
829 {
830 sprintf(buf,"%s","--");
831 }
832 else
833 {
834 /* Round to two decimals by printing. */
835 sprintf(tmp,"%.1e",ee);
836 sscanf(tmp,"%lf",&rnd);
837 sprintf(buf,"%g",rnd);
838 }
839
840 return buf;
841 }
842
843 static void analyse_ener(gmx_bool bCorr,const char *corrfn,
844 gmx_bool bFee,gmx_bool bSum,gmx_bool bFluct,gmx_bool bTempFluct,
845 gmx_bool bVisco,const char *visfn,int nmol,
846 int nconstr,
847 gmx_large_int_t start_step,double start_t,
848 gmx_large_int_t step,double t,
849 double time[], real reftemp,
850 enerdata_t *edat,
851 int nset,int set[],gmx_bool *bIsEner,
852 char **leg,gmx_enxnm_t *enm,
853 real Vaver,real ezero,
854 int nbmin,int nbmax,
855 const output_env_t oenv)
856 {
857 FILE *fp;
858 /* Check out the printed manual for equations! */
859 double Dt,aver,stddev,errest,delta_t,totaldrift;
860 enerdata_t *esum=NULL;
861 real xxx,integral,intBulk;
862 real sfrac,oldfrac,diffsum,diffav,fstep,pr_aver,pr_stddev,pr_errest;
863 double beta=0,expE,expEtot,*fee=NULL;
864 gmx_large_int_t nsteps;
865 int nexact,nnotexact;
866 double x1m,x1mk;
867 real Temp=-1,Pres=-1,VarV=-1,VarT=-1,VarEtot=-1,AvEtot=0,VarEnthalpy=-1;
868 int i,j,nout;
869 real chi2;
870 char buf[256],eebuf[100];
871
872 nsteps = step - start_step + 1;
873 if (nsteps < 1) {
874 fprintf(stdout,"Not enough steps (%s) for statistics\n",
875 gmx_step_str(nsteps,buf));
876 }
877 else {
878 /* Calculate the time difference */
879 delta_t = t - start_t;
880
881 fprintf(stdout,"\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
882 gmx_step_str(nsteps,buf),start_t,t,nset);
883
884 calc_averages(nset,edat,nbmin,nbmax);
885
886 if (bSum) {
887 esum = calc_sum(nset,edat,nbmin,nbmax);
888 }
889
890 if (edat->npoints == 0) {
891 nexact = 0;
892 nnotexact = nset;
893 } else {
894 nexact = 0;
895 nnotexact = 0;
896 for(i=0; (i<nset); i++) {
897 if (edat->s[i].bExactStat) {
898 nexact++;
899 } else {
900 nnotexact++;
901 }
902 }
903 }
904
905 if (nnotexact == 0) {
906 fprintf(stdout,"All statistics are over %s points\n",
907 gmx_step_str(edat->npoints,buf));
908 } else if (nexact == 0 || edat->npoints == edat->nframes) {
909 fprintf(stdout,"All statistics are over %d points (frames)\n",
910 edat->nframes);
911 } else {
912 fprintf(stdout,"The term%s",nnotexact==1 ? "" : "s");
913 for(i=0; (i<nset); i++) {
914 if (!edat->s[i].bExactStat) {
915 fprintf(stdout," '%s'",leg[i]);
916 }
917 }
918 fprintf(stdout," %s has statistics over %d points (frames)\n",
919 nnotexact==1 ? "is" : "are",edat->nframes);
920 fprintf(stdout,"All other statistics are over %s points\n",
921 gmx_step_str(edat->npoints,buf));
922 }
923 fprintf(stdout,"\n");
924
925 fprintf(stdout,"%-24s %10s %10s %10s %10s",
926 "Energy","Average","Err.Est.","RMSD","Tot-Drift");
927 if (bFee)
928 fprintf(stdout," %10s\n","-kT ln<e^(E/kT)>");
929 else
930 fprintf(stdout,"\n");
931 fprintf(stdout,"-------------------------------------------------------------------------------\n");
932
933 /* Initiate locals, only used with -sum */
934 expEtot=0;
935 if (bFee) {
936 beta = 1.0/(BOLTZ*reftemp);
937 snew(fee,nset);
938 }
939 for(i=0; (i<nset); i++) {
940 aver = edat->s[i].av;
941 stddev = edat->s[i].rmsd;
942 errest = edat->s[i].ee;
943
944 if (bFee) {
945 expE = 0;
946 for(j=0; (j<edat->nframes); j++) {
947 expE += exp(beta*(edat->s[i].ener[j] - aver)/nmol);
948 }
949 if (bSum)
950 expEtot+=expE/edat->nframes;
951
952 fee[i] = log(expE/edat->nframes)/beta + aver/nmol;
953 }
954 if (strstr(leg[i],"empera") != NULL) {
955 VarT = sqr(stddev);
956 Temp = aver;
957 } else if (strstr(leg[i],"olum") != NULL) {
958 VarV = sqr(stddev);
959 Vaver= aver;
960 } else if (strstr(leg[i],"essure") != NULL) {
961 Pres = aver;
962 } else if (strstr(leg[i],"otal") != NULL) {
963 VarEtot = sqr(stddev);
964 AvEtot = aver;
965 } else if (strstr(leg[i],"nthalpy") != NULL) {
966 VarEnthalpy = sqr(stddev);
967 }
968 if (bIsEner[i]) {
969 pr_aver = aver/nmol-ezero;
970 pr_stddev = stddev/nmol;
971 pr_errest = errest/nmol;
972 }
973 else {
974 pr_aver = aver;
975 pr_stddev = stddev;
976 pr_errest = errest;
977 }
978
979 /* Multiply the slope in steps with the number of steps taken */
980 totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
981 if (bIsEner[i])
982 {
983 totaldrift /= nmol;
984 }
985
986 fprintf(stdout,"%-24s %10g %10s %10g %10g",
987 leg[i],pr_aver,ee_pr(pr_errest,eebuf),pr_stddev,totaldrift);
988 if (bFee)
989 fprintf(stdout," %10g",fee[i]);
990
991 fprintf(stdout," (%s)\n",enm[set[i]].unit);
992
993 if (bFluct) {
994 for(j=0; (j<edat->nframes); j++)
995 edat->s[i].ener[j] -= aver;
996 }
997 }
998 if (bSum) {
999 totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
1000 fprintf(stdout,"%-24s %10g %10s %10s %10g (%s)",
1001 "Total",esum->s[0].av/nmol,ee_pr(esum->s[0].ee/nmol,eebuf),
1002 "--",totaldrift/nmol,enm[set[0]].unit);
1003 /* pr_aver,pr_stddev,a,totaldrift */
1004 if (bFee)
1005 fprintf(stdout," %10g %10g\n",
1006 log(expEtot)/beta + esum->s[0].av/nmol,log(expEtot)/beta);
1007 else
1008 fprintf(stdout,"\n");
1009 }
1010 if (bTempFluct && Temp != -1) {
1011 printf("\nTemperature dependent fluctuation properties at T = %g. #constr/mol = %d\n",Temp,nconstr);
1012 if (nmol < 2)
1013 printf("Warning: nmol = %d, this may not be what you want.\n",
1014 nmol);
1015 if (VarV != -1) {
1016 real tmp = VarV/(Vaver*BOLTZ*Temp*PRESFAC);
1017
1018 printf("Isothermal Compressibility: %10g /%s\n",
1019 tmp,unit_pres_bar);
1020 printf("Adiabatic bulk modulus: %10g %s\n",
1021 1.0/tmp,unit_pres_bar);
1022 }
1023 if (VarEnthalpy != -1) {
1024 real Cp = 1000*((VarEnthalpy/nmol)/(BOLTZ*Temp*Temp) -
1025 0.5*BOLTZ*nconstr);
1026 printf("Heat capacity at constant pressure Cp: %10g J/mol K\n",Cp);
1027 }
1028 if ((VarV != -1) && (VarEnthalpy != -1)) {
1029 real aP = (sqrt(VarEnthalpy*VarV/nmol))/(BOLTZ*Vaver*Temp*Temp);
1030 printf("Thermal expansion coefficient alphaP: %10g 1/K\n",aP);
1031 }
1032 if ((VarV == -1) && (VarEtot != -1)) {
1033 real Cv = 1000*((VarEtot/nmol)/(BOLTZ*Temp*Temp) -
1034 0.5*BOLTZ*nconstr);
1035 printf("Heat capacity at constant volume Cv: %10g J/mol K\n",Cv);
1036 }
1037 please_cite(stdout,"Allen1987a");
1038 }
1039 /* Do correlation function */
1040 if (edat->nframes > 1)
1041 {
1042 Dt = delta_t/(edat->nframes - 1);
1043 }
1044 else
1045 {
1046 Dt = 0;
1047 }
1048 if (bVisco) {
1049 const char* leg[] = { "Shear", "Bulk" };
1050 real factor;
1051 real **eneset;
1052 real **enesum;
1053
1054 /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1055
1056 /* Symmetrise tensor! (and store in first three elements)
1057 * And subtract average pressure!
1058 */
1059 snew(eneset,12);
1060 for(i=0; i<12; i++) {
1061 snew(eneset[i],edat->nframes);
1062 }
1063 snew(enesum,3);
1064 for(i=0; i<3; i++) {
1065 snew(enesum[i],edat->nframes);
1066 }
1067 for(i=0; (i<edat->nframes); i++) {
1068 eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
1069 eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
1070 eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
1071 for(j=3; j<=11; j++) {
1072 eneset[j][i] = edat->s[j].ener[i];
1073 }
1074 eneset[11][i] -= Pres;
1075 enesum[0][i] = 0.5*(edat->s[1].es[i].sum+edat->s[3].es[i].sum);
1076 enesum[1][i] = 0.5*(edat->s[2].es[i].sum+edat->s[6].es[i].sum);
1077 enesum[2][i] = 0.5*(edat->s[5].es[i].sum+edat->s[7].es[i].sum);
1078 }
1079
1080 einstein_visco("evisco.xvg","eviscoi.xvg",
1081 3,edat->nframes,enesum,Vaver,Temp,nsteps,time,oenv);
1082
1083 /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1084 /* Do it for shear viscosity */
1085 strcpy(buf,"Shear Viscosity");
1086 low_do_autocorr(corrfn,oenv,buf,edat->nframes,3,
1087 (edat->nframes+1)/2,eneset,Dt,
1088 eacNormal,1,TRUE,FALSE,FALSE,0.0,0.0,0,1);
1089
1090 /* Now for bulk viscosity */
1091 strcpy(buf,"Bulk Viscosity");
1092 low_do_autocorr(corrfn,oenv,buf,edat->nframes,1,
1093 (edat->nframes+1)/2,&(eneset[11]),Dt,
1094 eacNormal,1,TRUE,FALSE,FALSE,0.0,0.0,0,1);
1095
1096 factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
1097 fp=xvgropen(visfn,buf,"Time (ps)","\\8h\\4 (cp)",oenv);
1098 xvgr_legend(fp,asize(leg),leg,oenv);
1099
1100 /* Use trapezium rule for integration */
1101 integral = 0;
1102 intBulk = 0;
1103 nout = get_acfnout();
1104 if ((nout < 2) || (nout >= edat->nframes/2))
1105 nout = edat->nframes/2;
1106 for(i=1; (i<nout); i++)
1107 {
1108 integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor;
1109 intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1110 fprintf(fp,"%10g %10g %10g\n",(i*Dt),integral,intBulk);
1111 }
1112 ffclose(fp);
1113 }
1114 else if (bCorr) {
1115 if (bFluct)
1116 strcpy(buf,"Autocorrelation of Energy Fluctuations");
1117 else
1118 strcpy(buf,"Energy Autocorrelation");
1119 #if 0
1120 do_autocorr(corrfn,oenv,buf,edat->nframes,
1121 bSum ? 1 : nset,
1122 bSum ? &edat->s[nset-1].ener : eneset,
1123 (delta_t/edat->nframes),eacNormal,FALSE);
1124 #endif
1125 }
1126 }
1127 }
1128
1129 static void print_time(FILE *fp,double t)
1130 {
1131 fprintf(fp,"%12.6f",t);
1132 }
1133
1134 static void print1(FILE *fp,gmx_bool bDp,real e)
1135 {
1136 if (bDp)
1137 fprintf(fp," %16.12f",e);
1138 else
1139 fprintf(fp," %10.6f",e);
1140 }
1141
1142 static void fec(const char *ene2fn, const char *runavgfn,
1143 real reftemp, int nset, int set[], char *leg[],
1144 enerdata_t *edat, double time[],
1145 const output_env_t oenv)
1146 {
1147 const char* ravgleg[] = { "\\8D\\4E = E\\sB\\N-E\\sA\\N",
1148 "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N" };
1149 FILE *fp;
1150 ener_file_t enx;
1151 int nre,timecheck,step,nenergy,nenergy2,maxenergy;
1152 int i,j;
1153 gmx_bool bCont;
1154 real aver, beta;
1155 real **eneset2;
1156 double dE, sum;
1157 gmx_enxnm_t *enm=NULL;
1158 t_enxframe *fr;
1159 char buf[22];
1160
1161 /* read second energy file */
1162 snew(fr,1);
1163 enm = NULL;
1164 enx = open_enx(ene2fn,"r");
1165 do_enxnms(enx,&(fr->nre),&enm);
1166
1167 snew(eneset2,nset+1);
1168 nenergy2=0;
1169 maxenergy=0;
1170 timecheck=0;
1171 do {
1172 /* This loop searches for the first frame (when -b option is given),
1173 * or when this has been found it reads just one energy frame
1174 */
1175 do {
1176 bCont = do_enx(enx,fr);
1177
1178 if (bCont)
1179 timecheck = check_times(fr->t);
1180
1181 } while (bCont && (timecheck < 0));
1182
1183 /* Store energies for analysis afterwards... */
1184 if ((timecheck == 0) && bCont) {
1185 if (fr->nre > 0) {
1186 if ( nenergy2 >= maxenergy ) {
1187 maxenergy += 1000;
1188 for(i=0; i<=nset; i++)
1189 srenew(eneset2[i],maxenergy);
1190 }
1191 if (fr->t != time[nenergy2])
1192 fprintf(stderr,"\nWARNING time mismatch %g!=%g at frame %s\n",
1193 fr->t, time[nenergy2], gmx_step_str(fr->step,buf));
1194 for(i=0; i<nset; i++)
1195 eneset2[i][nenergy2] = fr->ener[set[i]].e;
1196 nenergy2++;
1197 }
1198 }
1199 } while (bCont && (timecheck == 0));
1200
1201 /* check */
1202 if (edat->nframes != nenergy2) {
1203 fprintf(stderr,"\nWARNING file length mismatch %d!=%d\n",
1204 edat->nframes,nenergy2);
1205 }
1206 nenergy = min(edat->nframes,nenergy2);
1207
1208 /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1209 fp=NULL;
1210 if (runavgfn) {
1211 fp=xvgropen(runavgfn,"Running average free energy difference",
1212 "Time (" unit_time ")","\\8D\\4E (" unit_energy ")",oenv);
1213 xvgr_legend(fp,asize(ravgleg),ravgleg,oenv);
1214 }
1215 fprintf(stdout,"\n%-24s %10s\n",
1216 "Energy","dF = -kT ln < exp(-(EB-EA)/kT) >A");
1217 sum=0;
1218 beta = 1.0/(BOLTZ*reftemp);
1219 for(i=0; i<nset; i++) {
1220 if (gmx_strcasecmp(leg[i],enm[set[i]].name)!=0)
1221 fprintf(stderr,"\nWARNING energy set name mismatch %s!=%s\n",
1222 leg[i],enm[set[i]].name);
1223 for(j=0; j<nenergy; j++) {
1224 dE = eneset2[i][j] - edat->s[i].ener[j];
1225 sum += exp(-dE*beta);
1226 if (fp)
1227 fprintf(fp,"%10g %10g %10g\n",
1228 time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
1229 }
1230 aver = -BOLTZ*reftemp*log(sum/nenergy);
1231 fprintf(stdout,"%-24s %10g\n",leg[i],aver);
1232 }
1233 if(fp) ffclose(fp);
1234 sfree(fr);
1235 }
1236
1237
1238 static void do_dhdl(t_enxframe *fr, FILE **fp_dhdl, const char *filename,
1239 int *blocks, int *hists, int *samples, int *nlambdas,
1240 const output_env_t oenv)
1241 {
1242 const char *dhdl="dH/d\\lambda",*deltag="\\DeltaH",*lambda="\\lambda";
1243 char title[STRLEN],label_x[STRLEN],label_y[STRLEN], legend[STRLEN];
1244 char buf[STRLEN];
1245 gmx_bool first=FALSE;
1246 int nblock_hist=0, nblock_dh=0, nblock_dhcoll=0;
1247 int i,j,k;
1248 /* coll data */
1249 double temp=0, start_time=0, delta_time=0, start_lambda=0, delta_lambda=0;
1250 static int setnr=0;
1251
1252 /* now count the blocks & handle the global dh data */
1253 for(i=0;i<fr->nblock;i++)
1254 {
1255 if (fr->block[i].id == enxDHHIST)
1256 {
1257 nblock_hist++;
1258 }
1259 else if (fr->block[i].id == enxDH)
1260 {
1261 nblock_dh++;
1262 }
1263 else if (fr->block[i].id == enxDHCOLL)
1264 {
1265 nblock_dhcoll++;
1266 if ( (fr->block[i].nsub < 1) ||
1267 (fr->block[i].sub[0].type != xdr_datatype_double) ||
1268 (fr->block[i].sub[0].nr < 5))
1269 {
1270 gmx_fatal(FARGS, "Unexpected block data");
1271 }
1272
1273 /* read the data from the DHCOLL block */
1274 temp = fr->block[i].sub[0].dval[0];
1275 start_time = fr->block[i].sub[0].dval[1];
1276 delta_time = fr->block[i].sub[0].dval[2];
1277 start_lambda = fr->block[i].sub[0].dval[3];
1278 delta_lambda = fr->block[i].sub[0].dval[4];
1279 }
1280 }
1281
1282 if (nblock_hist == 0 && nblock_dh == 0)
1283 {
1284 /* don't do anything */
1285 return;
1286 }
1287 if (nblock_hist>0 && nblock_dh>0)
1288 {
1289 gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
1290 }
1291 if (! *fp_dhdl )
1292 {
1293 if (nblock_dh>0)
1294 {
1295 sprintf(title,"%s, %s",dhdl,deltag);
1296 sprintf(label_x,"%s (%s)","Time",unit_time);
1297 sprintf(label_y,"(%s)",unit_energy);
1298 }
1299 else
1300 {
1301 sprintf(title,"N(%s)",deltag);
1302 sprintf(label_x,"%s (%s)",deltag,unit_energy);
1303 sprintf(label_y,"Samples");
1304 }
1305 *fp_dhdl=xvgropen_type(filename, title, label_x, label_y, exvggtXNY,
1306 oenv);
1307 sprintf(buf,"T = %g (K), %s = %g", temp, lambda, start_lambda);
1308 xvgr_subtitle(*fp_dhdl,buf,oenv);
1309 first=TRUE;
1310 }
1311
1312
1313
1314 (*hists)+=nblock_hist;
1315 (*blocks)+=nblock_dh;
1316 (*nlambdas) = nblock_hist+nblock_dh;
1317
1318
1319 /* write the data */
1320 if (nblock_hist > 0)
1321 {
1322 gmx_large_int_t sum=0;
1323 /* histograms */
1324 for(i=0;i<fr->nblock;i++)
1325 {
1326 t_enxblock *blk=&(fr->block[i]);
1327 if (blk->id==enxDHHIST)
1328 {
1329 double foreign_lambda, dx;
1330 gmx_large_int_t x0;
1331 int nhist, derivative;
1332
1333 /* check the block types etc. */
1334 if ( (blk->nsub < 2) ||
1335 (blk->sub[0].type != xdr_datatype_double) ||
1336 (blk->sub[1].type != xdr_datatype_large_int) ||
1337 (blk->sub[0].nr < 2) ||
1338 (blk->sub[1].nr < 2) )
1339 {
1340 gmx_fatal(FARGS, "Unexpected block data in file");
1341 }
1342 foreign_lambda=blk->sub[0].dval[0];
1343 dx=blk->sub[0].dval[1];
1344 nhist=blk->sub[1].lval[0];
1345 derivative=blk->sub[1].lval[1];
1346 for(j=0;j<nhist;j++)
1347 {
1348 const char *lg[1];
1349 x0=blk->sub[1].lval[2+j];
1350
1351 if (!derivative)
1352 {
1353 sprintf(legend,
1354 "N(%s(%s=%g) | %s=%g)",
1355 deltag, lambda, foreign_lambda,
1356 lambda, start_lambda);
1357 }
1358 else
1359 {
1360 sprintf(legend, "N(%s | %s=%g)",
1361 dhdl, lambda, start_lambda);
1362 }
1363
1364 lg[0]=legend;
1365 xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1366 setnr++;
1367 for(k=0;k<blk->sub[j+2].nr;k++)
1368 {
1369 int hist;
1370 double xmin, xmax;
1371
1372 hist=blk->sub[j+2].ival[k];
1373 xmin=(x0+k)*dx;
1374 xmax=(x0+k+1)*dx;
1375 fprintf(*fp_dhdl,"%g %d\n%g %d\n", xmin, hist,
1376 xmax, hist);
1377 sum+=hist;
1378 }
1379 /* multiple histogram data blocks in one histogram
1380 mean that the second one is the reverse of the first one:
1381 for dhdl derivatives, it's important to know both the
1382 maximum and minimum values */
1383 dx=-dx;
1384 }
1385 }
1386 }
1387
1388 (*samples) += (int)(sum/nblock_hist);
1389 }
1390 else
1391 {
1392 /* raw dh */
1393 int len=0;
1394 char **setnames=NULL;
1395 int nnames=nblock_dh;
1396
1397 if (fabs(delta_lambda) > 1e-9)
1398 {
1399 nnames++;
1400 }
1401 if (first)
1402 {
1403 snew(setnames, nnames);
1404 }
1405 j=0;
1406
1407 if (fabs(delta_lambda) > 1e-9)
1408 {
1409 /* lambda is a plotted value */
1410 setnames[j]=gmx_strdup(lambda);
1411 j++;
1412 }
1413
1414
1415 for(i=0;i<fr->nblock;i++)
1416 {
1417 t_enxblock *blk=&(fr->block[i]);
1418 if (blk->id == enxDH)
1419 {
1420 if (first)
1421 {
1422 /* do the legends */
1423 int derivative;
1424 double foreign_lambda;
1425
1426 derivative=blk->sub[0].ival[0];
1427 foreign_lambda=blk->sub[1].dval[0];
1428
1429 if (derivative)
1430 {
1431 sprintf(buf, "%s %s %g",dhdl,lambda,start_lambda);
1432 }
1433 else
1434 {
1435 sprintf(buf, "%s %s %g",deltag,lambda,
1436 foreign_lambda);
1437 }
1438 setnames[j] = gmx_strdup(buf);
1439 j++;
1440 }
1441
1442 if (len == 0)
1443 {
1444 len=blk->sub[2].nr;
1445 }
1446 else
1447 {
1448 if (len!=blk->sub[2].nr)
1449 {
1450 gmx_fatal(FARGS,
1451 "Length inconsistency in dhdl data");
1452 }
1453 }
1454 }
1455 }
1456
1457
1458 if (first)
1459 {
1460 xvgr_legend(*fp_dhdl, nblock_dh, (const char**)setnames, oenv);
1461 setnr += nblock_dh;
1462 for(i=0;i<nblock_dh;i++)
1463 {
1464 sfree(setnames[i]);
1465 }
1466 sfree(setnames);
1467 }
1468
1469 (*samples) += len;
1470 for(i=0;i<len;i++)
1471 {
1472 double time=start_time + delta_time*i;
1473
1474 fprintf(*fp_dhdl,"%.4f", time);
1475 if (fabs(delta_lambda) > 1e-9)
1476 {
1477 double lambda_now=i*delta_lambda + start_lambda;
1478 fprintf(*fp_dhdl," %.4f", lambda_now);
1479 }
1480 for(j=0;j<fr->nblock;j++)
1481 {
1482 t_enxblock *blk=&(fr->block[j]);
1483 if (blk->id == enxDH)
1484 {
1485 double value;
1486 if (blk->sub[2].type == xdr_datatype_float)
1487 {
1488 value=blk->sub[2].fval[i];
1489 }
1490 else
1491 {
1492 value=blk->sub[2].dval[i];
1493 }
1494 fprintf(*fp_dhdl," %g", value);
1495 }
1496 }
1497 fprintf(*fp_dhdl, "\n");
1498 }
1499 }
1500 }
1501
1502
1503 int gmx_energy(int argc,char *argv[])
1504 {
1505 const char *desc[] = {
1506
1507 "g_energy extracts energy components or distance restraint",
1508 "data from an energy file. The user is prompted to interactively",
1509 "select the energy terms she wants.[PAR]",
1510
1511 "Average, RMSD and drift are calculated with full precision from the",
1512 "simulation (see printed manual). Drift is calculated by performing",
1513 "a LSQ fit of the data to a straight line. The reported total drift",
1514 "is the difference of the fit at the first and last point.",
1515 "An error estimate of the average is given based on a block averages",
1516 "over 5 blocks using the full precision averages. The error estimate",
1517 "can be performed over multiple block lengths with the options",
1518 "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1519 "Note that in most cases the energy files contains averages over all",
1520 "MD steps, or over many more points than the number of frames in",
1521 "energy file. This makes the g_energy statistics output more accurate",
1522 "than the xvg output. When exact averages are not present in the energy",
1523 "file the statistics mentioned above is simply over the single, per-frame",
1524 "energy values.[PAR]",
1525
1526 "The term fluctuation gives the RMSD around the LSQ fit.[PAR]",
1527
1528 "Some fluctuation-dependent properties can be calculated provided",
1529 "the correct energy terms are selected. The following properties",
1530 "will be computed:[BR]",
1531 "Property Energy terms needed[BR]",
1532 "---------------------------------------------------[BR]",
1533 "Heat capacity Cp (NPT sims): Enthalpy, Temp [BR]",
1534 "Heat capacity Cv (NVT sims): Etot, Temp [BR]",
1535 "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp[BR]",
1536 "Isothermal compressibility: Vol, Temp [BR]",
1537 "Adiabatic bulk modulus: Vol, Temp [BR]",
1538 "---------------------------------------------------[BR]",
1539 "You always need to set the number of molecules [TT]-nmol[tt], and,",
1540 "if you used constraints in your simulations you will need to give",
1541 "the number of constraints per molecule [TT]-nconstr[tt] in order to",
1542 "correct for this: (nconstr/2) kB is subtracted from the heat",
1543 "capacity in this case. For instance in the case of rigid water",
1544 "you need to give the value 3 to this option.[PAR]",
1545
1546 "When the [TT]-viol[tt] option is set, the time averaged",
1547 "violations are plotted and the running time-averaged and",
1548 "instantaneous sum of violations are recalculated. Additionally",
1549 "running time-averaged and instantaneous distances between",
1550 "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1551
1552 "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1553 "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1554 "The first two options plot the orientation, the last three the",
1555 "deviations of the orientations from the experimental values.",
1556 "The options that end on an 'a' plot the average over time",
1557 "as a function of restraint. The options that end on a 't'",
1558 "prompt the user for restraint label numbers and plot the data",
1559 "as a function of time. Option [TT]-odr[tt] plots the RMS",
1560 "deviation as a function of restraint.",
1561 "When the run used time or ensemble averaged orientation restraints,",
1562 "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1563 "not ensemble-averaged orientations and deviations instead of",
1564 "the time and ensemble averages.[PAR]",
1565
1566 "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1567 "tensor for each orientation restraint experiment. With option",
1568 "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1569
1570 "Option [TT]-odh[tt] extracts and plots the free energy data",
1571 "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1572 "from the ener.edr file.[PAR]",
1573
1574 "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1575 "difference with an ideal gas state: [BR]",
1576 " Delta A = A(N,V,T) - A_idgas(N,V,T) = kT ln < e^(Upot/kT) >[BR]",
1577 " Delta G = G(N,p,T) - G_idgas(N,p,T) = kT ln < e^(Upot/kT) >[BR]",
1578 "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1579 "the average is over the ensemble (or time in a trajectory).",
1580 "Note that this is in principle",
1581 "only correct when averaging over the whole (Boltzmann) ensemble",
1582 "and using the potential energy. This also allows for an entropy",
1583 "estimate using:[BR]",
1584 " Delta S(N,V,T) = S(N,V,T) - S_idgas(N,V,T) = (<Upot> - Delta A)/T[BR]",
1585 " Delta S(N,p,T) = S(N,p,T) - S_idgas(N,p,T) = (<Upot> + pV - Delta G)/T",
1586 "[PAR]",
1587
1588 "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1589 "difference is calculated dF = -kT ln < e ^ -(EB-EA)/kT >A ,",
1590 "where EA and EB are the energies from the first and second energy",
1591 "files, and the average is over the ensemble A. [BB]NOTE[bb] that",
1592 "the energies must both be calculated from the same trajectory."
1593
1594 };
1595 static gmx_bool bSum=FALSE,bFee=FALSE,bPrAll=FALSE,bFluct=FALSE;
1596 static gmx_bool bDp=FALSE,bMutot=FALSE,bOrinst=FALSE,bOvec=FALSE;
1597 static int skip=0,nmol=1,nconstr=0,nbmin=5,nbmax=5;
1598 static real reftemp=300.0,ezero=0;
1599 t_pargs pa[] = {
1600 { "-fee", FALSE, etBOOL, {&bFee},
1601 "Do a free energy estimate" },
1602 { "-fetemp", FALSE, etREAL,{&reftemp},
1603 "Reference temperature for free energy calculation" },
1604 { "-zero", FALSE, etREAL, {&ezero},
1605 "Subtract a zero-point energy" },
1606 { "-sum", FALSE, etBOOL, {&bSum},
1607 "Sum the energy terms selected rather than display them all" },
1608 { "-dp", FALSE, etBOOL, {&bDp},
1609 "Print energies in high precision" },
1610 { "-nbmin", FALSE, etINT, {&nbmin},
1611 "Minimum number of blocks for error estimate" },
1612 { "-nbmax", FALSE, etINT, {&nbmax},
1613 "Maximum number of blocks for error estimate" },
1614 { "-mutot",FALSE, etBOOL, {&bMutot},
1615 "Compute the total dipole moment from the components" },
1616 { "-skip", FALSE, etINT, {&skip},
1617 "Skip number of frames between data points" },
1618 { "-aver", FALSE, etBOOL, {&bPrAll},
1619 "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1620 { "-nmol", FALSE, etINT, {&nmol},
1621 "Number of molecules in your sample: the energies are divided by this number" },
1622 { "-nconstr", FALSE, etINT, {&nconstr},
1623 "Number of constraints per molecule. Necessary for calculating the heat capacity" },
1624 { "-fluc", FALSE, etBOOL, {&bFluct},
1625 "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1626 { "-orinst", FALSE, etBOOL, {&bOrinst},
1627 "Analyse instantaneous orientation data" },
1628 { "-ovec", FALSE, etBOOL, {&bOvec},
1629 "Also plot the eigenvectors with -oten" }
1630 };
1631 const char* drleg[] = {
1632 "Running average",
1633 "Instantaneous"
1634 };
1635 static const char *setnm[] = {
1636 "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1637 "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1638 "Volume", "Pressure"
1639 };
1640
1641 FILE *out=NULL,*fp_pairs=NULL,*fort=NULL,*fodt=NULL,*foten=NULL;
1642 FILE *fp_dhdl=NULL;
1643 FILE **drout;
1644 ener_file_t fp;
1645 int timecheck=0;
1646 gmx_mtop_t mtop;
1647 gmx_localtop_t *top=NULL;
1648 t_inputrec ir;
1649 t_energy **ee;
1650 enerdata_t edat;
1651 gmx_enxnm_t *enm=NULL;
1652 t_enxframe *frame,*fr=NULL;
1653 int cur=0;
1654 #define NEXT (1-cur)
1655 int nre,teller,teller_disre,nfr;
1656 gmx_large_int_t start_step;
1657 int nor=0,nex=0,norfr=0,enx_i=0;
1658 real start_t;
1659 real *bounds=NULL,*violaver=NULL,*oobs=NULL,*orient=NULL,*odrms=NULL;
1660 int *index=NULL,*pair=NULL,norsel=0,*orsel=NULL,*or_label=NULL;
1661 int nbounds=0,npairs;
1662 gmx_bool bDisRe,bDRAll,bORA,bORT,bODA,bODR,bODT,bORIRE,bOTEN,bDHDL;
1663 gmx_bool bFoundStart,bCont,bEDR,bVisco;
1664 double sum,sumaver,sumt,ener,dbl;
1665 double *time=NULL;
1666 real Vaver;
1667 int *set=NULL,i,j,k,nset,sss;
1668 gmx_bool *bIsEner=NULL;
1669 char **pairleg,**odtleg,**otenleg;
1670 char **leg=NULL;
1671 char **nms;
1672 char *anm_j,*anm_k,*resnm_j,*resnm_k;
1673 int resnr_j,resnr_k;
1674 const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
1675 char buf[256];
1676 output_env_t oenv;
1677 t_enxblock *blk=NULL;
1678 t_enxblock *blk_disre=NULL;
1679 int ndisre=0;
1680 int dh_blocks=0, dh_hists=0, dh_samples=0, dh_lambdas=0;
1681
1682 t_filenm fnm[] = {
1683 { efEDR, "-f", NULL, ffREAD },
1684 { efEDR, "-f2", NULL, ffOPTRD },
1685 { efTPX, "-s", NULL, ffOPTRD },
1686 { efXVG, "-o", "energy", ffWRITE },
1687 { efXVG, "-viol", "violaver",ffOPTWR },
1688 { efXVG, "-pairs","pairs", ffOPTWR },
1689 { efXVG, "-ora", "orienta", ffOPTWR },
1690 { efXVG, "-ort", "orientt", ffOPTWR },
1691 { efXVG, "-oda", "orideva", ffOPTWR },
1692 { efXVG, "-odr", "oridevr", ffOPTWR },
1693 { efXVG, "-odt", "oridevt", ffOPTWR },
1694 { efXVG, "-oten", "oriten", ffOPTWR },
1695 { efXVG, "-corr", "enecorr", ffOPTWR },
1696 { efXVG, "-vis", "visco", ffOPTWR },
1697 { efXVG, "-ravg", "runavgdf",ffOPTWR },
1698 { efXVG, "-odh", "dhdl" ,ffOPTWR }
1699 };
1700 #define NFILE asize(fnm)
1701 int npargs;
1702 t_pargs *ppa;
1703
1704 CopyRight(stderr,argv[0]);
1705 npargs = asize(pa);
1706 ppa = add_acf_pargs(&npargs,pa);
1707 parse_common_args(&argc,argv,
1708 PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
1709 NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
1710
1711 bDRAll = opt2bSet("-pairs",NFILE,fnm);
1712 bDisRe = opt2bSet("-viol",NFILE,fnm) || bDRAll;
1713 bORA = opt2bSet("-ora",NFILE,fnm);
1714 bORT = opt2bSet("-ort",NFILE,fnm);
1715 bODA = opt2bSet("-oda",NFILE,fnm);
1716 bODR = opt2bSet("-odr",NFILE,fnm);
1717 bODT = opt2bSet("-odt",NFILE,fnm);
1718 bORIRE = bORA || bORT || bODA || bODR || bODT;
1719 bOTEN = opt2bSet("-oten",NFILE,fnm);
1720 bDHDL = opt2bSet("-odh",NFILE,fnm);
1721
1722 nset = 0;
1723
1724 snew(frame,2);
1725 fp = open_enx(ftp2fn(efEDR,NFILE,fnm),"r");
1726 do_enxnms(fp,&nre,&enm);
1727
1728 Vaver = -1;
1729
1730 bVisco = opt2bSet("-vis",NFILE,fnm);
1731
1732 if (!bDisRe && !bDHDL)
1733 {
1734 if (bVisco) {
1735 nset=asize(setnm);
1736 snew(set,nset);
1737 /* This is nasty code... To extract Pres tensor, Volume and Temperature */
1738 for(j=0; j<nset; j++) {
1739 for(i=0; i<nre; i++) {
1740 if (strstr(enm[i].name,setnm[j])) {
1741 set[j]=i;
1742 break;
1743 }
1744 }
1745 if (i == nre) {
1746 if (gmx_strcasecmp(setnm[j],"Volume")==0) {
1747 printf("Enter the box volume (" unit_volume "): ");
1748 if(1 != scanf("%lf",&dbl))
1749 {
1750 gmx_fatal(FARGS,"Error reading user input");
1751 }
1752 Vaver = dbl;
1753 } else
1754 gmx_fatal(FARGS,"Could not find term %s for viscosity calculation",
1755 setnm[j]);
1756 }
1757 }
1758 }
1759 else
1760 {
1761 set=select_by_name(nre,enm,&nset);
1762 }
1763 /* Print all the different units once */
1764 sprintf(buf,"(%s)",enm[set[0]].unit);
1765 for(i=1; i<nset; i++) {
1766 for(j=0; j<i; j++) {
1767 if (strcmp(enm[set[i]].unit,enm[set[j]].unit) == 0) {
1768 break;
1769 }
1770 }
1771 if (j == i) {
1772 strcat(buf,", (");
1773 strcat(buf,enm[set[i]].unit);
1774 strcat(buf,")");
1775 }
1776 }
1777 out=xvgropen(opt2fn("-o",NFILE,fnm),"Gromacs Energies","Time (ps)",buf,
1778 oenv);
1779
1780 snew(leg,nset+1);
1781 for(i=0; (i<nset); i++)
1782 leg[i] = enm[set[i]].name;
1783 if (bSum) {
1784 leg[nset]=strdup("Sum");
1785 xvgr_legend(out,nset+1,(const char**)leg,oenv);
1786 }
1787 else
1788 xvgr_legend(out,nset,(const char**)leg,oenv);
1789
1790 snew(bIsEner,nset);
1791 for(i=0; (i<nset); i++) {
1792 bIsEner[i] = FALSE;
1793 for (j=0; (j <= F_ETOT); j++)
1794 bIsEner[i] = bIsEner[i] ||
1795 (gmx_strcasecmp(interaction_function[j].longname,leg[i]) == 0);
1796 }
1797
1798 if (bPrAll && nset > 1) {
1799 gmx_fatal(FARGS,"Printing averages can only be done when a single set is selected");
1800 }
1801
1802 time = NULL;
1803
1804 if (bORIRE || bOTEN)
1805 get_orires_parms(ftp2fn(efTPX,NFILE,fnm),&nor,&nex,&or_label,&oobs);
1806
1807 if (bORIRE) {
1808 if (bOrinst)
1809 enx_i = enxORI;
1810 else
1811 enx_i = enxOR;
1812
1813 if (bORA || bODA)
1814 snew(orient,nor);
1815 if (bODR)
1816 snew(odrms,nor);
1817 if (bORT || bODT) {
1818 fprintf(stderr,"Select the orientation restraint labels you want (-1 is all)\n");
1819 fprintf(stderr,"End your selection with 0\n");
1820 j = -1;
1821 orsel = NULL;
1822 do {
1823 j++;
1824 srenew(orsel,j+1);
1825 if(1 != scanf("%d",&(orsel[j])))
1826 {
1827 gmx_fatal(FARGS,"Error reading user input");
1828 }
1829 } while (orsel[j] > 0);
1830 if (orsel[0] == -1) {
1831 fprintf(stderr,"Selecting all %d orientation restraints\n",nor);
1832 norsel = nor;
1833 srenew(orsel,nor);
1834 for(i=0; i<nor; i++)
1835 orsel[i] = i;
1836 } else {
1837 /* Build the selection */
1838 norsel=0;
1839 for(i=0; i<j; i++) {
1840 for(k=0; k<nor; k++)
1841 if (or_label[k] == orsel[i]) {
1842 orsel[norsel] = k;
1843 norsel++;
1844 break;
1845 }
1846 if (k == nor)
1847 fprintf(stderr,"Orientation restraint label %d not found\n",
1848 orsel[i]);
1849 }
1850 }
1851 snew(odtleg,norsel);
1852 for(i=0; i<norsel; i++) {
1853 snew(odtleg[i],256);
1854 sprintf(odtleg[i],"%d",or_label[orsel[i]]);
1855 }
1856 if (bORT) {
1857 fort=xvgropen(opt2fn("-ort",NFILE,fnm), "Calculated orientations",
1858 "Time (ps)","",oenv);
1859 if (bOrinst)
1860 fprintf(fort,"%s",orinst_sub);
1861 xvgr_legend(fort,norsel,(const char**)odtleg,oenv);
1862 }
1863 if (bODT) {
1864 fodt=xvgropen(opt2fn("-odt",NFILE,fnm),
1865 "Orientation restraint deviation",
1866 "Time (ps)","",oenv);
1867 if (bOrinst)
1868 fprintf(fodt,"%s",orinst_sub);
1869 xvgr_legend(fodt,norsel,(const char**)odtleg,oenv);
1870 }
1871 }
1872 }
1873 if (bOTEN) {
1874 foten=xvgropen(opt2fn("-oten",NFILE,fnm),
1875 "Order tensor","Time (ps)","",oenv);
1876 snew(otenleg,bOvec ? nex*12 : nex*3);
1877 for(i=0; i<nex; i++) {
1878 for(j=0; j<3; j++) {
1879 sprintf(buf,"eig%d",j+1);
1880 otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
1881 }
1882 if (bOvec) {
1883 for(j=0; j<9; j++) {
1884 sprintf(buf,"vec%d%s",j/3+1,j%3==0 ? "x" : (j%3==1 ? "y" : "z"));
1885 otenleg[12*i+3+j] = strdup(buf);
1886 }
1887 }
1888 }
1889 xvgr_legend(foten,bOvec ? nex*12 : nex*3,(const char**)otenleg,oenv);
1890 }
1891 }
1892 else if (bDisRe)
1893 {
1894 nbounds=get_bounds(ftp2fn(efTPX,NFILE,fnm),&bounds,&index,&pair,&npairs,
1895 &mtop,&top,&ir);
1896 snew(violaver,npairs);
1897 out=xvgropen(opt2fn("-o",NFILE,fnm),"Sum of Violations",
1898 "Time (ps)","nm",oenv);
1899 xvgr_legend(out,2,drleg,oenv);
1900 if (bDRAll) {
1901 fp_pairs=xvgropen(opt2fn("-pairs",NFILE,fnm),"Pair Distances",
1902 "Time (ps)","Distance (nm)",oenv);
1903 if (output_env_get_print_xvgr_codes(oenv))
1904 fprintf(fp_pairs,"@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
1905 ir.dr_tau);
1906 }
1907 }
1908
1909
1910 /* Initiate energies and set them to zero */
1911 edat.nsteps = 0;
1912 edat.npoints = 0;
1913 edat.nframes = 0;
1914 edat.step = NULL;
1915 edat.steps = NULL;
1916 edat.points = NULL;
1917 snew(edat.s,nset);
1918
1919 /* Initiate counters */
1920 teller = 0;
1921 teller_disre = 0;
1922 bFoundStart = FALSE;
1923 start_step = 0;
1924 start_t = 0;
1925 do {
1926 /* This loop searches for the first frame (when -b option is given),
1927 * or when this has been found it reads just one energy frame
1928 */
1929 do {
1930 bCont = do_enx(fp,&(frame[NEXT]));
1931
1932 if (bCont) {
1933 timecheck = check_times(frame[NEXT].t);
1934 }
1935 } while (bCont && (timecheck < 0));
1936
1937 if ((timecheck == 0) && bCont) {
1938 /* We read a valid frame, so we can use it */
1939 fr = &(frame[NEXT]);
1940
1941 if (fr->nre > 0) {
1942 /* The frame contains energies, so update cur */
1943 cur = NEXT;
1944
1945 if (edat.nframes % 1000 == 0)
1946 {
1947 srenew(edat.step,edat.nframes+1000);
1948 srenew(edat.steps,edat.nframes+1000);
1949 srenew(edat.points,edat.nframes+1000);
1950 for(i=0; i<nset; i++)
1951 {
1952 srenew(edat.s[i].ener,edat.nframes+1000);
1953 srenew(edat.s[i].es ,edat.nframes+1000);
1954 }
1955 }
1956
1957 nfr = edat.nframes;
1958 edat.step[nfr] = fr->step;
1959
1960 if (!bFoundStart)
1961 {
1962 bFoundStart = TRUE;
1963 /* Initiate the previous step data */
1964 start_step = fr->step;
1965 start_t = fr->t;
1966 /* Initiate the energy sums */
1967 edat.steps[nfr] = 1;
1968 edat.points[nfr] = 1;
1969 for(i=0; i<nset; i++)
1970 {
1971 sss = set[i];
1972 edat.s[i].es[nfr].sum = fr->ener[sss].e;
1973 edat.s[i].es[nfr].sum2 = 0;
1974 }
1975 edat.nsteps = 1;
1976 edat.npoints = 1;
1977 }
1978 else
1979 {
1980 edat.steps[nfr] = fr->nsteps;
1981 {
1982 if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
1983 {
1984 if (fr->nsum <= 1)
1985 {
1986 edat.points[nfr] = 1;
1987 for(i=0; i<nset; i++)
1988 {
1989 sss = set[i];
1990 edat.s[i].es[nfr].sum = fr->ener[sss].e;
1991 edat.s[i].es[nfr].sum2 = 0;
1992 }
1993 edat.npoints += 1;
1994 }
1995 else
1996 {
1997 edat.points[nfr] = fr->nsum;
1998 for(i=0; i<nset; i++)
1999 {
2000 sss = set[i];
2001 edat.s[i].es[nfr].sum = fr->ener[sss].esum;
2002 edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
2003 }
2004 edat.npoints += fr->nsum;
2005 }
2006 }
2007 else
2008 {
2009 /* The interval does not match fr->nsteps:
2010 * can not do exact averages.
2011 */
2012 edat.npoints = 0;
2013 }
2014 edat.nsteps = fr->step - start_step + 1;
2015 }
2016 }
2017 for(i=0; i<nset; i++)
2018 {
2019 edat.s[i].ener[nfr] = fr->ener[set[i]].e;
2020 }
2021 }
2022 /*
2023 * Define distance restraint legends. Can only be done after
2024 * the first frame has been read... (Then we know how many there are)
2025 */
2026 blk_disre=find_block_id_enxframe(fr, enxDISRE, NULL);
2027 if (bDisRe && bDRAll && !leg && blk_disre)
2028 {
2029 t_iatom *fa;
2030 t_iparams *ip;
2031
2032 fa = top->idef.il[F_DISRES].iatoms;
2033 ip = top->idef.iparams;
2034 if (blk_disre->nsub != 2 ||
2035 (blk_disre->sub[0].nr != blk_disre->sub[1].nr) )
2036 {
2037 gmx_incons("Number of disre sub-blocks not equal to 2");
2038 }
2039
2040 ndisre=blk_disre->sub[0].nr ;
2041 if (ndisre != top->idef.il[F_DISRES].nr/3)
2042 {
2043 gmx_fatal(FARGS,"Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
2044 ndisre,top->idef.il[F_DISRES].nr/3);
2045 }
2046 snew(pairleg,ndisre);
2047 for(i=0; i<ndisre; i++)
2048 {
2049 snew(pairleg[i],30);
2050 j=fa[3*i+1];
2051 k=fa[3*i+2];
2052 gmx_mtop_atominfo_global(&mtop,j,&anm_j,&resnr_j,&resnm_j);
2053 gmx_mtop_atominfo_global(&mtop,k,&anm_k,&resnr_k,&resnm_k);
2054 sprintf(pairleg[i],"%d %s %d %s (%d)",
2055 resnr_j,anm_j,resnr_k,anm_k,
2056 ip[fa[3*i]].disres.label);
2057 }
2058 set=select_it(ndisre,pairleg,&nset);
2059 snew(leg,2*nset);
2060 for(i=0; (i<nset); i++)
2061 {
2062 snew(leg[2*i],32);
2063 sprintf(leg[2*i], "a %s",pairleg[set[i]]);
2064 snew(leg[2*i+1],32);
2065 sprintf(leg[2*i+1],"i %s",pairleg[set[i]]);
2066 }
2067 xvgr_legend(fp_pairs,2*nset,(const char**)leg,oenv);
2068 }
2069
2070 /*
2071 * Store energies for analysis afterwards...
2072 */
2073 if (!bDisRe && !bDHDL && (fr->nre > 0)) {
2074 if (edat.nframes % 1000 == 0) {
2075 srenew(time,edat.nframes+1000);
2076 }
2077 time[edat.nframes] = fr->t;
2078 edat.nframes++;
2079 }
2080 /*
2081 * Printing time, only when we do not want to skip frames
2082 */
2083 if (!skip || teller % skip == 0) {
2084 if (bDisRe) {
2085 /*******************************************
2086 * D I S T A N C E R E S T R A I N T S
2087 *******************************************/
2088 if (ndisre > 0)
2089 {
2090 #ifndef GMX_DOUBLE
2091 float *disre_rt = blk_disre->sub[0].fval;
2092 float *disre_rm3tav = blk_disre->sub[1].fval;
2093 #else
2094 double *disre_rt = blk_disre->sub[0].dval;
2095 double *disre_rm3tav = blk_disre->sub[1].dval;
2096 #endif
2097
2098 print_time(out,fr->t);
2099 if (violaver == NULL)
2100 snew(violaver,ndisre);
2101
2102 /* Subtract bounds from distances, to calculate violations */
2103 calc_violations(disre_rt, disre_rm3tav,
2104 nbounds,pair,bounds,violaver,&sumt,&sumaver);
2105
2106 fprintf(out," %8.4f %8.4f\n",sumaver,sumt);
2107 if (bDRAll) {
2108 print_time(fp_pairs,fr->t);
2109 for(i=0; (i<nset); i++) {
2110 sss=set[i];
2111 fprintf(fp_pairs," %8.4f", mypow(disre_rm3tav[sss],minthird));
2112 fprintf(fp_pairs," %8.4f", disre_rt[sss]);
2113 }
2114 fprintf(fp_pairs,"\n");
2115 }
2116 teller_disre++;
2117 }
2118 }
2119 else if (bDHDL)
2120 {
2121 do_dhdl(fr, &fp_dhdl, opt2fn("-odh",NFILE,fnm),
2122 &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas,
2123 oenv);
2124 }
2125 /*******************************************
2126 * E N E R G I E S
2127 *******************************************/
2128 else {
2129 if (fr->nre > 0) {
2130 if (bPrAll)
2131 {
2132 /* We skip frames with single points (usually only the first frame),
2133 * since they would result in an average plot with outliers.
2134 */
2135 if (fr->nsum > 1) {
2136 print_time(out,fr->t);
2137 print1(out,bDp,fr->ener[set[0]].e);
2138 print1(out,bDp,fr->ener[set[0]].esum/fr->nsum);
2139 print1(out,bDp,sqrt(fr->ener[set[0]].eav/fr->nsum));
2140 fprintf(out,"\n");
2141 }
2142 }
2143 else
2144 {
2145 print_time(out,fr->t);
2146 if (bSum)
2147 {
2148 sum = 0;
2149 for(i=0; i<nset; i++)
2150 {
2151 sum += fr->ener[set[i]].e;
2152 }
2153 print1(out,bDp,sum/nmol-ezero);
2154 }
2155 else
2156 {
2157 for(i=0; (i<nset); i++)
2158 {
2159 if (bIsEner[i])
2160 {
2161 print1(out,bDp,(fr->ener[set[i]].e)/nmol-ezero);
2162 }
2163 else
2164 {
2165 print1(out,bDp,fr->ener[set[i]].e);
2166 }
2167 }
2168 }
2169 fprintf(out,"\n");
2170 }
2171 }
2172 #if 0
2173 /* we first count the blocks that have id 0: the orire blocks */
2174 block_orire=0;
2175 for(b=0;b<fr->nblock;b++)
2176 {
2177 if (fr->block[b].id == mde_block_type_orire)
2178 nblock_orire++;
2179 }
2180 #endif
2181 blk = find_block_id_enxframe(fr, enx_i, NULL);
2182 if (bORIRE && blk)
2183 {
2184 #ifndef GMX_DOUBLE
2185 xdr_datatype dt=xdr_datatype_float;
2186 #else
2187 xdr_datatype dt=xdr_datatype_double;
2188 #endif
2189 real *vals;
2190
2191 if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
2192 {
2193 gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
2194 }
2195 #ifndef GMX_DOUBLE
2196 vals=blk->sub[0].fval;
2197 #else
2198 vals=blk->sub[0].dval;
2199 #endif
2200
2201 if (blk->sub[0].nr != (size_t)nor)
2202 gmx_fatal(FARGS,"Number of orientation restraints in energy file (%d) does not match with the topology (%d)", blk->sub[0].nr);
2203 if (bORA || bODA)
2204 {
2205 for(i=0; i<nor; i++)
2206 orient[i] += vals[i];
2207 }
2208 if (bODR)
2209 {
2210 for(i=0; i<nor; i++)
2211 odrms[i] += sqr(vals[i]-oobs[i]);
2212 }
2213 if (bORT)
2214 {
2215 fprintf(fort," %10f",fr->t);
2216 for(i=0; i<norsel; i++)
2217 fprintf(fort," %g",vals[orsel[i]]);
2218 fprintf(fort,"\n");
2219 }
2220 if (bODT)
2221 {
2222 fprintf(fodt," %10f",fr->t);
2223 for(i=0; i<norsel; i++)
2224 fprintf(fodt," %g", vals[orsel[i]]-oobs[orsel[i]]);
2225 fprintf(fodt,"\n");
2226 }
2227 norfr++;
2228 }
2229 blk = find_block_id_enxframe(fr, enxORT, NULL);
2230 if (bOTEN && blk)
2231 {
2232 #ifndef GMX_DOUBLE
2233 xdr_datatype dt=xdr_datatype_float;
2234 #else
2235 xdr_datatype dt=xdr_datatype_double;
2236 #endif
2237 real *vals;
2238
2239 if ( (blk->nsub != 1) || (blk->sub[0].type!=dt) )
2240 gmx_fatal(FARGS,"Orientational restraints read in incorrectly");
2241 #ifndef GMX_DOUBLE
2242 vals=blk->sub[0].fval;
2243 #else
2244 vals=blk->sub[0].dval;
2245 #endif
2246
2247 if (blk->sub[0].nr != (size_t)(nex*12))
2248 gmx_fatal(FARGS,"Number of orientation experiments in energy file (%g) does not match with the topology (%d)",
2249 blk->sub[0].nr/12, nex);
2250 fprintf(foten," %10f",fr->t);
2251 for(i=0; i<nex; i++)
2252 for(j=0; j<(bOvec?12:3); j++)
2253 fprintf(foten," %g",vals[i*12+j]);
2254 fprintf(foten,"\n");
2255 }
2256 }
2257 }
2258 }
2259 } while (bCont && (timecheck == 0));
2260
2261 fprintf(stderr,"\n");
2262 close_enx(fp);
2263 if (out)
2264 ffclose(out);
2265
2266 if (bDRAll)
2267 ffclose(fp_pairs);
2268
2269 if (bORT)
2270 ffclose(fort);
2271 if (bODT)
2272 ffclose(fodt);
2273 if (bORA)
2274 {
2275 out = xvgropen(opt2fn("-ora",NFILE,fnm),
2276 "Average calculated orientations",
2277 "Restraint label","",oenv);
2278 if (bOrinst)
2279 fprintf(out,"%s",orinst_sub);
2280 for(i=0; i<nor; i++)
2281 fprintf(out,"%5d %g\n",or_label[i],orient[i]/norfr);
2282 ffclose(out);
2283 }
2284 if (bODA) {
2285 out = xvgropen(opt2fn("-oda",NFILE,fnm),
2286 "Average restraint deviation",
2287 "Restraint label","",oenv);
2288 if (bOrinst)
2289 fprintf(out,"%s",orinst_sub);
2290 for(i=0; i<nor; i++)
2291 fprintf(out,"%5d %g\n",or_label[i],orient[i]/norfr-oobs[i]);
2292 ffclose(out);
2293 }
2294 if (bODR) {
2295 out = xvgropen(opt2fn("-odr",NFILE,fnm),
2296 "RMS orientation restraint deviations",
2297 "Restraint label","",oenv);
2298 if (bOrinst)
2299 fprintf(out,"%s",orinst_sub);
2300 for(i=0; i<nor; i++)
2301 fprintf(out,"%5d %g\n",or_label[i],sqrt(odrms[i]/norfr));
2302 ffclose(out);
2303 }
2304 if (bOTEN)
2305 ffclose(foten);
2306
2307 if (bDisRe)
2308 {
2309 analyse_disre(opt2fn("-viol",NFILE,fnm),
2310 teller_disre,violaver,bounds,index,pair,nbounds,oenv);
2311 }
2312 else if (bDHDL)
2313 {
2314 if (fp_dhdl)
2315 {
2316 ffclose(fp_dhdl);
2317 printf("\n\nWrote %d lambda values with %d samples as ",
2318 dh_lambdas, dh_samples);
2319 if (dh_hists > 0)
2320 {
2321 printf("%d dH histograms ", dh_hists);
2322 }
2323 if (dh_blocks> 0)
2324 {
2325 printf("%d dH data blocks ", dh_blocks);
2326 }
2327 printf("to %s\n", opt2fn("-odh",NFILE,fnm));
2328
2329 }
2330 else
2331 {
2332 gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f",NFILE,fnm));
2333 }
2334
2335 }
2336 else
2337 {
2338 analyse_ener(opt2bSet("-corr",NFILE,fnm),opt2fn("-corr",NFILE,fnm),
2339 bFee,bSum,bFluct,opt2parg_bSet("-nmol",npargs,ppa),
2340 bVisco,opt2fn("-vis",NFILE,fnm),
2341 nmol,nconstr,start_step,start_t,frame[cur].step,frame[cur].t,
2342 time,reftemp,&edat,
2343 nset,set,bIsEner,leg,enm,Vaver,ezero,nbmin,nbmax,
2344 oenv);
2345 }
2346 if (opt2bSet("-f2",NFILE,fnm)) {
2347 fec(opt2fn("-f2",NFILE,fnm), opt2fn("-ravg",NFILE,fnm),
2348 reftemp, nset, set, leg, &edat, time ,oenv);
2349 }
2350
2351 {
2352 const char *nxy = "-nxy";
2353
2354 do_view(oenv,opt2fn("-o",NFILE,fnm),nxy);
2355 do_view(oenv,opt2fn_null("-ravg",NFILE,fnm),nxy);
2356 do_view(oenv,opt2fn_null("-ora",NFILE,fnm),nxy);
2357 do_view(oenv,opt2fn_null("-ort",NFILE,fnm),nxy);
2358 do_view(oenv,opt2fn_null("-oda",NFILE,fnm),nxy);
2359 do_view(oenv,opt2fn_null("-odr",NFILE,fnm),nxy);
2360 do_view(oenv,opt2fn_null("-odt",NFILE,fnm),nxy);
2361 do_view(oenv,opt2fn_null("-oten",NFILE,fnm),nxy);
2362 do_view(oenv,opt2fn_null("-odh",NFILE,fnm),nxy);
2363 }
2364 thanx(stderr);
2365
2366 return 0;
2367 }
got merge a long time ago