?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Fix for MSVC linking and release-build OpenMM plugin dir.
[gromacs/qmmm-gamess-us.git] / src / tools / gmx_energy.c
blob552090887c6817a261f7f3933b5c2292cf921d35
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 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 bool *bE;
101 int n,k,j,i;
102 int *set;
103 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 bool *bE;
160 int n,k,kk,j,i,nmatch,nind,nss;
161 int *set;
162 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 (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 (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 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(bool bCorr,const char *corrfn,
844 bool bFee,bool bSum,bool bFluct,bool bTempFluct,
845 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[],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;
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 thru %.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 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 for(i=1; (i<edat->nframes/2); i++) {
1104 integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor;
1105 intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1106 fprintf(fp,"%10g %10g %10g\n",(i*Dt),integral,intBulk);
1107 }
1108 ffclose(fp);
1109 }
1110 else if (bCorr) {
1111 if (bFluct)
1112 strcpy(buf,"Autocorrelation of Energy Fluctuations");
1113 else
1114 strcpy(buf,"Energy Autocorrelation");
1115 #if 0
1116 do_autocorr(corrfn,oenv,buf,edat->nframes,
1117 bSum ? 1 : nset,
1118 bSum ? &edat->s[nset-1].ener : eneset,
1119 (delta_t/edat->nframes),eacNormal,FALSE);
1120 #endif
1121 }
1122 }
1123 }
1124
1125 static void print_time(FILE *fp,double t)
1126 {
1127 fprintf(fp,"%12.6f",t);
1128 }
1129
1130 static void print1(FILE *fp,bool bDp,real e)
1131 {
1132 if (bDp)
1133 fprintf(fp," %16.12f",e);
1134 else
1135 fprintf(fp," %10.6f",e);
1136 }
1137
1138 static void fec(const char *ene2fn, const char *runavgfn,
1139 real reftemp, int nset, int set[], char *leg[],
1140 enerdata_t *edat, double time[],
1141 const output_env_t oenv)
1142 {
1143 char *ravgleg[] = { "\\8D\\4E = E\\sB\\N-E\\sA\\N",
1144 "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N" };
1145 FILE *fp;
1146 ener_file_t enx;
1147 int nre,timecheck,step,nenergy,nenergy2,maxenergy;
1148 int i,j;
1149 bool bCont;
1150 real aver, beta;
1151 real **eneset2;
1152 double dE, sum;
1153 gmx_enxnm_t *enm=NULL;
1154 t_enxframe *fr;
1155 char buf[22];
1156
1157 /* read second energy file */
1158 snew(fr,1);
1159 enm = NULL;
1160 enx = open_enx(ene2fn,"r");
1161 do_enxnms(enx,&(fr->nre),&enm);
1162
1163 snew(eneset2,nset+1);
1164 nenergy2=0;
1165 maxenergy=0;
1166 timecheck=0;
1167 do {
1168 /* This loop searches for the first frame (when -b option is given),
1169 * or when this has been found it reads just one energy frame
1170 */
1171 do {
1172 bCont = do_enx(enx,fr);
1173
1174 if (bCont)
1175 timecheck = check_times(fr->t);
1176
1177 } while (bCont && (timecheck < 0));
1178
1179 /* Store energies for analysis afterwards... */
1180 if ((timecheck == 0) && bCont) {
1181 if (fr->nre > 0) {
1182 if ( nenergy2 >= maxenergy ) {
1183 maxenergy += 1000;
1184 for(i=0; i<=nset; i++)
1185 srenew(eneset2[i],maxenergy);
1186 }
1187 if (fr->t != time[nenergy2])
1188 fprintf(stderr,"\nWARNING time mismatch %g!=%g at frame %s\n",
1189 fr->t, time[nenergy2], gmx_step_str(fr->step,buf));
1190 for(i=0; i<nset; i++)
1191 eneset2[i][nenergy2] = fr->ener[set[i]].e;
1192 nenergy2++;
1193 }
1194 }
1195 } while (bCont && (timecheck == 0));
1196
1197 /* check */
1198 if (edat->nframes != nenergy2) {
1199 fprintf(stderr,"\nWARNING file length mismatch %d!=%d\n",
1200 edat->nframes,nenergy2);
1201 }
1202 nenergy = min(edat->nframes,nenergy2);
1203
1204 /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1205 fp=NULL;
1206 if (runavgfn) {
1207 fp=xvgropen(runavgfn,"Running average free energy difference",
1208 "Time (" unit_time ")","\\8D\\4E (" unit_energy ")",oenv);
1209 xvgr_legend(fp,asize(ravgleg),ravgleg,oenv);
1210 }
1211 fprintf(stdout,"\n%-24s %10s\n",
1212 "Energy","dF = -kT ln < exp(-(EB-EA)/kT) >A");
1213 sum=0;
1214 beta = 1.0/(BOLTZ*reftemp);
1215 for(i=0; i<nset; i++) {
1216 if (strcasecmp(leg[i],enm[set[i]].name)!=0)
1217 fprintf(stderr,"\nWARNING energy set name mismatch %s!=%s\n",
1218 leg[i],enm[set[i]].name);
1219 for(j=0; j<nenergy; j++) {
1220 dE = eneset2[i][j] - edat->s[i].ener[j];
1221 sum += exp(-dE*beta);
1222 if (fp)
1223 fprintf(fp,"%10g %10g %10g\n",
1224 time[j], dE, -BOLTZ*reftemp*log(sum/(j+1)) );
1225 }
1226 aver = -BOLTZ*reftemp*log(sum/nenergy);
1227 fprintf(stdout,"%-24s %10g\n",leg[i],aver);
1228 }
1229 if(fp) ffclose(fp);
1230 sfree(fr);
1231 }
1232
1233 int gmx_energy(int argc,char *argv[])
1234 {
1235 const char *desc[] = {
1236
1237 "g_energy extracts energy components or distance restraint",
1238 "data from an energy file. The user is prompted to interactively",
1239 "select the energy terms she wants.[PAR]",
1240
1241 "Average, RMSD and drift are calculated with full precision from the",
1242 "simulation (see printed manual). Drift is calculated by performing",
1243 "a LSQ fit of the data to a straight line. The reported total drift",
1244 "is the difference of the fit at the first and last point.",
1245 "An error estimate of the average is given based on a block averages",
1246 "over 5 blocks using the full precision averages. The error estimate",
1247 "can be performed over multiple block lengths with the options",
1248 "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1249 "Note that in most cases the energy files contains averages over all",
1250 "MD steps, or over many more points than the number of frames in",
1251 "energy file. This makes the g_energy statistics output more accurate",
1252 "than the xvg output. When exact averages are not present in the energy",
1253 "file the statistics mentioned above is simply over the single, per-frame",
1254 "energy values.[PAR]",
1255
1256 "The term fluctuation gives the RMSD around the LSQ fit.[PAR]",
1257
1258 "Some fluctuation-dependent properties can be calculated provided",
1259 "the correct energy terms are selected. The following properties",
1260 "will be computed:[BR]",
1261 "Property Energy terms needed[BR]",
1262 "--------------------------------------------------[BR]",
1263 "Heat capacity Cp (NPT sims) Enthalpy, Temp [BR]",
1264 "Heat capacity Cv (NVT sims) Etot, Temp [BR]",
1265 "Thermal expansion coeff. (NPT) Enthalpy, Vol, Temp[BR]",
1266 "Isothermal compressibility Vol, Temp [BR]",
1267 "Adiabatic bulk modulus Vol, Temp [PBR]",
1268 "--------------------------------------------------[BR]",
1269 "You always need to set the number of molecules [TT]-nmol[tt], and,",
1270 "if you used constraints in your simulations you will need to give",
1271 "the number of constraints per molecule [TT]-nconstr[tt] in order to",
1272 "correct for this: (nconstr/2) kB is subtracted from the heat",
1273 "capacity in this case. For instance in the case of rigid water",
1274 "you need to give the value 3 to this option.[PAR]",
1275
1276 "When the [TT]-viol[tt] option is set, the time averaged",
1277 "violations are plotted and the running time-averaged and",
1278 "instantaneous sum of violations are recalculated. Additionally",
1279 "running time-averaged and instantaneous distances between",
1280 "selected pairs can be plotted with the [TT]-pairs[tt] option.[PAR]",
1281
1282 "Options [TT]-ora[tt], [TT]-ort[tt], [TT]-oda[tt], [TT]-odr[tt] and",
1283 "[TT]-odt[tt] are used for analyzing orientation restraint data.",
1284 "The first two options plot the orientation, the last three the",
1285 "deviations of the orientations from the experimental values.",
1286 "The options that end on an 'a' plot the average over time",
1287 "as a function of restraint. The options that end on a 't'",
1288 "prompt the user for restraint label numbers and plot the data",
1289 "as a function of time. Option [TT]-odr[tt] plots the RMS",
1290 "deviation as a function of restraint.",
1291 "When the run used time or ensemble averaged orientation restraints,",
1292 "option [TT]-orinst[tt] can be used to analyse the instantaneous,",
1293 "not ensemble-averaged orientations and deviations instead of",
1294 "the time and ensemble averages.[PAR]",
1295
1296 "Option [TT]-oten[tt] plots the eigenvalues of the molecular order",
1297 "tensor for each orientation restraint experiment. With option",
1298 "[TT]-ovec[tt] also the eigenvectors are plotted.[PAR]",
1299
1300 "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1301 "difference with an ideal gas state: [BR]",
1302 " Delta A = A(N,V,T) - A_idgas(N,V,T) = kT ln < e^(Upot/kT) >[BR]",
1303 " Delta G = G(N,p,T) - G_idgas(N,p,T) = kT ln < e^(Upot/kT) >[BR]",
1304 "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and"
1305 "the average is over the ensemble (or time in a trajectory).",
1306 "Note that this is in principle",
1307 "only correct when averaging over the whole (Boltzmann) ensemble",
1308 "and using the potential energy. This also allows for an entropy",
1309 "estimate using:[BR]",
1310 " Delta S(N,V,T) = S(N,V,T) - S_idgas(N,V,T) = (<Upot> - Delta A)/T[BR]",
1311 " Delta S(N,p,T) = S(N,p,T) - S_idgas(N,p,T) = (<Upot> + pV - Delta G)/T",
1312 "[PAR]",
1313
1314 "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1315 "difference is calculated dF = -kT ln < e ^ -(EB-EA)/kT >A ,",
1316 "where EA and EB are the energies from the first and second energy",
1317 "files, and the average is over the ensemble A. [BB]NOTE[bb] that",
1318 "the energies must both be calculated from the same trajectory."
1319
1320 };
1321 static bool bSum=FALSE,bFee=FALSE,bPrAll=FALSE,bFluct=FALSE;
1322 static bool bDp=FALSE,bMutot=FALSE,bOrinst=FALSE,bOvec=FALSE;
1323 static int skip=0,nmol=1,nconstr=0,nbmin=5,nbmax=5;
1324 static real reftemp=300.0,ezero=0;
1325 t_pargs pa[] = {
1326 { "-fee", FALSE, etBOOL, {&bFee},
1327 "Do a free energy estimate" },
1328 { "-fetemp", FALSE, etREAL,{&reftemp},
1329 "Reference temperature for free energy calculation" },
1330 { "-zero", FALSE, etREAL, {&ezero},
1331 "Subtract a zero-point energy" },
1332 { "-sum", FALSE, etBOOL, {&bSum},
1333 "Sum the energy terms selected rather than display them all" },
1334 { "-dp", FALSE, etBOOL, {&bDp},
1335 "Print energies in high precision" },
1336 { "-nbmin", FALSE, etINT, {&nbmin},
1337 "Minimum number of blocks for error estimate" },
1338 { "-nbmax", FALSE, etINT, {&nbmax},
1339 "Maximum number of blocks for error estimate" },
1340 { "-mutot",FALSE, etBOOL, {&bMutot},
1341 "Compute the total dipole moment from the components" },
1342 { "-skip", FALSE, etINT, {&skip},
1343 "Skip number of frames between data points" },
1344 { "-aver", FALSE, etBOOL, {&bPrAll},
1345 "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1346 { "-nmol", FALSE, etINT, {&nmol},
1347 "Number of molecules in your sample: the energies are divided by this number" },
1348 { "-nconstr", FALSE, etINT, {&nconstr},
1349 "Number of constraints per molecule. Necessary for calculating the heat capacity" },
1350 { "-fluc", FALSE, etBOOL, {&bFluct},
1351 "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1352 { "-orinst", FALSE, etBOOL, {&bOrinst},
1353 "Analyse instantaneous orientation data" },
1354 { "-ovec", FALSE, etBOOL, {&bOvec},
1355 "Also plot the eigenvectors with -oten" }
1356 };
1357 char *drleg[] = {
1358 "Running average",
1359 "Instantaneous"
1360 };
1361 static const char *setnm[] = {
1362 "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1363 "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1364 "Volume", "Pressure"
1365 };
1366
1367 FILE *out,*fp_pairs=NULL,*fort=NULL,*fodt=NULL,*foten=NULL;
1368 FILE **drout;
1369 ener_file_t fp;
1370 int timecheck=0;
1371 gmx_mtop_t mtop;
1372 gmx_localtop_t *top=NULL;
1373 t_inputrec ir;
1374 t_energy **ee;
1375 enerdata_t edat;
1376 gmx_enxnm_t *enm=NULL;
1377 t_enxframe *frame,*fr=NULL;
1378 int cur=0;
1379 #define NEXT (1-cur)
1380 int nre,teller,teller_disre,nfr;
1381 gmx_large_int_t start_step;
1382 int nor=0,nex=0,norfr=0,enx_i=0;
1383 real start_t;
1384 real *bounds=NULL,*violaver=NULL,*oobs=NULL,*orient=NULL,*odrms=NULL;
1385 int *index=NULL,*pair=NULL,norsel=0,*orsel=NULL,*or_label=NULL;
1386 int nbounds=0,npairs;
1387 bool bDisRe,bDRAll,bORA,bORT,bODA,bODR,bODT,bORIRE,bOTEN;
1388 bool bFoundStart,bCont,bEDR,bVisco;
1389 double sum,sumaver,sumt,ener,dbl;
1390 double *time=NULL;
1391 real Vaver;
1392 int *set=NULL,i,j,k,nset,sss;
1393 bool *bIsEner=NULL;
1394 char **pairleg,**odtleg,**otenleg;
1395 char **leg=NULL;
1396 char **nms;
1397 char *anm_j,*anm_k,*resnm_j,*resnm_k;
1398 int resnr_j,resnr_k;
1399 const char *orinst_sub = "@ subtitle \"instantaneous\"\n";
1400 char buf[256];
1401 output_env_t oenv;
1402 t_filenm fnm[] = {
1403 { efEDR, "-f", NULL, ffREAD },
1404 { efEDR, "-f2", NULL, ffOPTRD },
1405 { efTPX, "-s", NULL, ffOPTRD },
1406 { efXVG, "-o", "energy", ffWRITE },
1407 { efXVG, "-viol", "violaver",ffOPTWR },
1408 { efXVG, "-pairs","pairs", ffOPTWR },
1409 { efXVG, "-ora", "orienta", ffOPTWR },
1410 { efXVG, "-ort", "orientt", ffOPTWR },
1411 { efXVG, "-oda", "orideva", ffOPTWR },
1412 { efXVG, "-odr", "oridevr", ffOPTWR },
1413 { efXVG, "-odt", "oridevt", ffOPTWR },
1414 { efXVG, "-oten", "oriten", ffOPTWR },
1415 { efXVG, "-corr", "enecorr", ffOPTWR },
1416 { efXVG, "-vis", "visco", ffOPTWR },
1417 { efXVG, "-ravg", "runavgdf",ffOPTWR }
1418 };
1419 #define NFILE asize(fnm)
1420 int npargs;
1421 t_pargs *ppa;
1422
1423 CopyRight(stderr,argv[0]);
1424 npargs = asize(pa);
1425 ppa = add_acf_pargs(&npargs,pa);
1426 parse_common_args(&argc,argv,
1427 PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END | PCA_BE_NICE,
1428 NFILE,fnm,npargs,ppa,asize(desc),desc,0,NULL,&oenv);
1429
1430 bDRAll = opt2bSet("-pairs",NFILE,fnm);
1431 bDisRe = opt2bSet("-viol",NFILE,fnm) || bDRAll;
1432 bORA = opt2bSet("-ora",NFILE,fnm);
1433 bORT = opt2bSet("-ort",NFILE,fnm);
1434 bODA = opt2bSet("-oda",NFILE,fnm);
1435 bODR = opt2bSet("-odr",NFILE,fnm);
1436 bODT = opt2bSet("-odt",NFILE,fnm);
1437 bORIRE = bORA || bORT || bODA || bODR || bODT;
1438 bOTEN = opt2bSet("-oten",NFILE,fnm);
1439
1440 nset = 0;
1441
1442 snew(frame,2);
1443 fp = open_enx(ftp2fn(efEDR,NFILE,fnm),"r");
1444 do_enxnms(fp,&nre,&enm);
1445
1446 Vaver = -1;
1447
1448 bVisco = opt2bSet("-vis",NFILE,fnm);
1449
1450 if (!bDisRe) {
1451 if (bVisco) {
1452 nset=asize(setnm);
1453 snew(set,nset);
1454 /* This is nasty code... To extract Pres tensor, Volume and Temperature */
1455 for(j=0; j<nset; j++) {
1456 for(i=0; i<nre; i++) {
1457 if (strstr(enm[i].name,setnm[j])) {
1458 set[j]=i;
1459 break;
1460 }
1461 }
1462 if (i == nre) {
1463 if (strcasecmp(setnm[j],"Volume")==0) {
1464 printf("Enter the box volume (" unit_volume "): ");
1465 if(1 != scanf("%lf",&dbl))
1466 {
1467 gmx_fatal(FARGS,"Error reading user input");
1468 }
1469 Vaver = dbl;
1470 } else
1471 gmx_fatal(FARGS,"Could not find term %s for viscosity calculation",
1472 setnm[j]);
1473 }
1474 }
1475 }
1476 else {
1477 set=select_by_name(nre,enm,&nset);
1478 }
1479 /* Print all the different units once */
1480 sprintf(buf,"(%s)",enm[set[0]].unit);
1481 for(i=1; i<nset; i++) {
1482 for(j=0; j<i; j++) {
1483 if (strcmp(enm[set[i]].unit,enm[set[j]].unit) == 0) {
1484 break;
1485 }
1486 }
1487 if (j == i) {
1488 strcat(buf,", (");
1489 strcat(buf,enm[set[i]].unit);
1490 strcat(buf,")");
1491 }
1492 }
1493 out=xvgropen(opt2fn("-o",NFILE,fnm),"Gromacs Energies","Time (ps)",buf,
1494 oenv);
1495
1496 snew(leg,nset+1);
1497 for(i=0; (i<nset); i++)
1498 leg[i] = enm[set[i]].name;
1499 if (bSum) {
1500 leg[nset]=strdup("Sum");
1501 xvgr_legend(out,nset+1,leg,oenv);
1502 }
1503 else
1504 xvgr_legend(out,nset,leg,oenv);
1505
1506 snew(bIsEner,nset);
1507 for(i=0; (i<nset); i++) {
1508 bIsEner[i] = FALSE;
1509 for (j=0; (j <= F_ETOT); j++)
1510 bIsEner[i] = bIsEner[i] ||
1511 (strcasecmp(interaction_function[j].longname,leg[i]) == 0);
1512 }
1513
1514 if (bPrAll && nset > 1) {
1515 gmx_fatal(FARGS,"Printing averages can only be done when a single set is selected");
1516 }
1517
1518 time = NULL;
1519
1520 if (bORIRE || bOTEN)
1521 get_orires_parms(ftp2fn(efTPX,NFILE,fnm),&nor,&nex,&or_label,&oobs);
1522
1523 if (bORIRE) {
1524 if (bOrinst)
1525 enx_i = enxORI;
1526 else
1527 enx_i = enxOR;
1528
1529 if (bORA || bODA)
1530 snew(orient,nor);
1531 if (bODR)
1532 snew(odrms,nor);
1533 if (bORT || bODT) {
1534 fprintf(stderr,"Select the orientation restraint labels you want (-1 is all)\n");
1535 fprintf(stderr,"End your selection with 0\n");
1536 j = -1;
1537 orsel = NULL;
1538 do {
1539 j++;
1540 srenew(orsel,j+1);
1541 if(1 != scanf("%d",&(orsel[j])))
1542 {
1543 gmx_fatal(FARGS,"Error reading user input");
1544 }
1545 } while (orsel[j] > 0);
1546 if (orsel[0] == -1) {
1547 fprintf(stderr,"Selecting all %d orientation restraints\n",nor);
1548 norsel = nor;
1549 srenew(orsel,nor);
1550 for(i=0; i<nor; i++)
1551 orsel[i] = i;
1552 } else {
1553 /* Build the selection */
1554 norsel=0;
1555 for(i=0; i<j; i++) {
1556 for(k=0; k<nor; k++)
1557 if (or_label[k] == orsel[i]) {
1558 orsel[norsel] = k;
1559 norsel++;
1560 break;
1561 }
1562 if (k == nor)
1563 fprintf(stderr,"Orientation restraint label %d not found\n",
1564 orsel[i]);
1565 }
1566 }
1567 snew(odtleg,norsel);
1568 for(i=0; i<norsel; i++) {
1569 snew(odtleg[i],256);
1570 sprintf(odtleg[i],"%d",or_label[orsel[i]]);
1571 }
1572 if (bORT) {
1573 fort=xvgropen(opt2fn("-ort",NFILE,fnm), "Calculated orientations",
1574 "Time (ps)","",oenv);
1575 if (bOrinst)
1576 fprintf(fort,"%s",orinst_sub);
1577 xvgr_legend(fort,norsel,odtleg,oenv);
1578 }
1579 if (bODT) {
1580 fodt=xvgropen(opt2fn("-odt",NFILE,fnm),
1581 "Orientation restraint deviation",
1582 "Time (ps)","",oenv);
1583 if (bOrinst)
1584 fprintf(fodt,"%s",orinst_sub);
1585 xvgr_legend(fodt,norsel,odtleg,oenv);
1586 }
1587 }
1588 }
1589 if (bOTEN) {
1590 foten=xvgropen(opt2fn("-oten",NFILE,fnm),
1591 "Order tensor","Time (ps)","",oenv);
1592 snew(otenleg,bOvec ? nex*12 : nex*3);
1593 for(i=0; i<nex; i++) {
1594 for(j=0; j<3; j++) {
1595 sprintf(buf,"eig%d",j+1);
1596 otenleg[(bOvec ? 12 : 3)*i+j] = strdup(buf);
1597 }
1598 if (bOvec) {
1599 for(j=0; j<9; j++) {
1600 sprintf(buf,"vec%d%s",j/3+1,j%3==0 ? "x" : (j%3==1 ? "y" : "z"));
1601 otenleg[12*i+3+j] = strdup(buf);
1602 }
1603 }
1604 }
1605 xvgr_legend(foten,bOvec ? nex*12 : nex*3,otenleg,oenv);
1606 }
1607 }
1608 else {
1609 nbounds=get_bounds(ftp2fn(efTPX,NFILE,fnm),&bounds,&index,&pair,&npairs,
1610 &mtop,&top,&ir);
1611 snew(violaver,npairs);
1612 out=xvgropen(opt2fn("-o",NFILE,fnm),"Sum of Violations",
1613 "Time (ps)","nm",oenv);
1614 xvgr_legend(out,2,drleg,oenv);
1615 if (bDRAll) {
1616 fp_pairs=xvgropen(opt2fn("-pairs",NFILE,fnm),"Pair Distances",
1617 "Time (ps)","Distance (nm)",oenv);
1618 if (output_env_get_print_xvgr_codes(oenv))
1619 fprintf(fp_pairs,"@ subtitle \"averaged (tau=%g) and instantaneous\"\n",
1620 ir.dr_tau);
1621 }
1622 }
1623
1624 /* Initiate energies and set them to zero */
1625 edat.nsteps = 0;
1626 edat.npoints = 0;
1627 edat.nframes = 0;
1628 edat.step = NULL;
1629 edat.steps = NULL;
1630 edat.points = NULL;
1631 snew(edat.s,nset);
1632
1633 /* Initiate counters */
1634 teller = 0;
1635 teller_disre = 0;
1636 bFoundStart = FALSE;
1637 start_step = 0;
1638 start_t = 0;
1639 do {
1640 /* This loop searches for the first frame (when -b option is given),
1641 * or when this has been found it reads just one energy frame
1642 */
1643 do {
1644 bCont = do_enx(fp,&(frame[NEXT]));
1645
1646 if (bCont) {
1647 timecheck = check_times(frame[NEXT].t);
1648 }
1649 } while (bCont && (timecheck < 0));
1650
1651 if ((timecheck == 0) && bCont) {
1652 /* We read a valid frame, so we can use it */
1653 fr = &(frame[NEXT]);
1654
1655 if (fr->nre > 0) {
1656 /* The frame contains energies, so update cur */
1657 cur = NEXT;
1658
1659 if (edat.nframes % 1000 == 0)
1660 {
1661 srenew(edat.step,edat.nframes+1000);
1662 srenew(edat.steps,edat.nframes+1000);
1663 srenew(edat.points,edat.nframes+1000);
1664 for(i=0; i<nset; i++)
1665 {
1666 srenew(edat.s[i].ener,edat.nframes+1000);
1667 srenew(edat.s[i].es ,edat.nframes+1000);
1668 }
1669 }
1670
1671 nfr = edat.nframes;
1672 edat.step[nfr] = fr->step;
1673
1674 if (!bFoundStart)
1675 {
1676 bFoundStart = TRUE;
1677 /* Initiate the previous step data */
1678 start_step = fr->step;
1679 start_t = fr->t;
1680 /* Initiate the energy sums */
1681 edat.steps[nfr] = 1;
1682 edat.points[nfr] = 1;
1683 for(i=0; i<nset; i++)
1684 {
1685 sss = set[i];
1686 edat.s[i].es[nfr].sum = fr->ener[sss].e;
1687 edat.s[i].es[nfr].sum2 = 0;
1688 }
1689 edat.nsteps = 1;
1690 edat.npoints = 1;
1691 }
1692 else
1693 {
1694 edat.steps[nfr] = fr->nsteps;
1695 {
1696 if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
1697 {
1698 if (fr->nsum <= 1)
1699 {
1700 edat.points[nfr] = 1;
1701 for(i=0; i<nset; i++)
1702 {
1703 sss = set[i];
1704 edat.s[i].es[nfr].sum = fr->ener[sss].e;
1705 edat.s[i].es[nfr].sum2 = 0;
1706 }
1707 edat.npoints += 1;
1708 }
1709 else
1710 {
1711 edat.points[nfr] = fr->nsum;
1712 for(i=0; i<nset; i++)
1713 {
1714 sss = set[i];
1715 edat.s[i].es[nfr].sum = fr->ener[sss].esum;
1716 edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
1717 }
1718 edat.npoints += fr->nsum;
1719 }
1720 }
1721 else
1722 {
1723 /* The interval does not match fr->nsteps:
1724 * can not do exact averages.
1725 */
1726 edat.npoints = 0;
1727 }
1728 edat.nsteps = fr->step - start_step + 1;
1729 }
1730 }
1731 for(i=0; i<nset; i++)
1732 {
1733 edat.s[i].ener[nfr] = fr->ener[set[i]].e;
1734 }
1735 }
1736 /*
1737 * Define distance restraint legends. Can only be done after
1738 * the first frame has been read... (Then we know how many there are)
1739 */
1740 if (bDisRe && bDRAll && !leg && (fr->ndisre > 0)) {
1741 t_iatom *fa;
1742 t_iparams *ip;
1743
1744 fa = top->idef.il[F_DISRES].iatoms;
1745 ip = top->idef.iparams;
1746
1747 if (fr->ndisre != top->idef.il[F_DISRES].nr/3)
1748 gmx_fatal(FARGS,"Number of disre pairs in the energy file (%d) does not match the number in the run input file (%d)\n",
1749 fr->ndisre,top->idef.il[F_DISRES].nr/3);
1750
1751 snew(pairleg,fr->ndisre);
1752 for(i=0; i<fr->ndisre; i++) {
1753 snew(pairleg[i],30);
1754 j=fa[3*i+1];
1755 k=fa[3*i+2];
1756 gmx_mtop_atominfo_global(&mtop,j,&anm_j,&resnr_j,&resnm_j);
1757 gmx_mtop_atominfo_global(&mtop,k,&anm_k,&resnr_k,&resnm_k);
1758 sprintf(pairleg[i],"%d %s %d %s (%d)",
1759 resnr_j+1,anm_j,resnr_k+1,anm_k,
1760 ip[fa[3*i]].disres.label);
1761 }
1762 set=select_it(fr->ndisre,pairleg,&nset);
1763 snew(leg,2*nset);
1764 for(i=0; (i<nset); i++) {
1765 snew(leg[2*i],32);
1766 sprintf(leg[2*i], "a %s",pairleg[set[i]]);
1767 snew(leg[2*i+1],32);
1768 sprintf(leg[2*i+1],"i %s",pairleg[set[i]]);
1769 }
1770 xvgr_legend(fp_pairs,2*nset,leg,oenv);
1771 }
1772
1773 /*
1774 * Store energies for analysis afterwards...
1775 */
1776 if (!bDisRe && (fr->nre > 0)) {
1777 if (edat.nframes % 1000 == 0) {
1778 srenew(time,edat.nframes+1000);
1779 }
1780 time[edat.nframes] = fr->t;
1781 edat.nframes++;
1782 }
1783 /*
1784 * Printing time, only when we do not want to skip frames
1785 */
1786 if (!skip || teller % skip == 0) {
1787 if (bDisRe) {
1788 /*******************************************
1789 * D I S T A N C E R E S T R A I N T S
1790 *******************************************/
1791 if (fr->ndisre > 0) {
1792 print_time(out,fr->t);
1793 if (violaver == NULL)
1794 snew(violaver,fr->ndisre);
1795
1796 /* Subtract bounds from distances, to calculate violations */
1797 calc_violations(fr->disre_rt,fr->disre_rm3tav,
1798 nbounds,pair,bounds,violaver,&sumt,&sumaver);
1799
1800 fprintf(out," %8.4f %8.4f\n",sumaver,sumt);
1801 if (bDRAll) {
1802 print_time(fp_pairs,fr->t);
1803 for(i=0; (i<nset); i++) {
1804 sss=set[i];
1805 fprintf(fp_pairs," %8.4f",
1806 mypow(fr->disre_rm3tav[sss],minthird));
1807 fprintf(fp_pairs," %8.4f",
1808 fr->disre_rt[sss]);
1809 }
1810 fprintf(fp_pairs,"\n");
1811 }
1812 teller_disre++;
1813 }
1814 }
1815 /*******************************************
1816 * E N E R G I E S
1817 *******************************************/
1818 else {
1819 if (fr->nre > 0) {
1820 if (bPrAll)
1821 {
1822 /* We skip frames with single points (usually only the first frame),
1823 * since they would result in an average plot with outliers.
1824 */
1825 if (fr->nsum > 1) {
1826 print_time(out,fr->t);
1827 print1(out,bDp,fr->ener[set[0]].e);
1828 print1(out,bDp,fr->ener[set[0]].esum/fr->nsum);
1829 print1(out,bDp,sqrt(fr->ener[set[0]].eav/fr->nsum));
1830 fprintf(out,"\n");
1831 }
1832 }
1833 else
1834 {
1835 print_time(out,fr->t);
1836 if (bSum)
1837 {
1838 sum = 0;
1839 for(i=0; i<nset; i++)
1840 {
1841 sum += fr->ener[set[i]].e;
1842 }
1843 print1(out,bDp,sum/nmol-ezero);
1844 }
1845 else
1846 {
1847 for(i=0; (i<nset); i++)
1848 {
1849 if (bIsEner[i])
1850 {
1851 print1(out,bDp,(fr->ener[set[i]].e)/nmol-ezero);
1852 }
1853 else
1854 {
1855 print1(out,bDp,fr->ener[set[i]].e);
1856 }
1857 }
1858 }
1859 fprintf(out,"\n");
1860 }
1861 }
1862 if (bORIRE && fr->nblock>enx_i && fr->nr[enx_i]>0) {
1863 if (fr->nr[enx_i] != nor)
1864 gmx_fatal(FARGS,"Number of orientation restraints in energy file (%d) does not match with the topology (%d)",fr->nr[enx_i],nor);
1865 if (bORA || bODA)
1866 for(i=0; i<nor; i++)
1867 orient[i] += fr->block[enx_i][i];
1868 if (bODR)
1869 for(i=0; i<nor; i++)
1870 odrms[i] += sqr(fr->block[enx_i][i]-oobs[i]);
1871 if (bORT) {
1872 fprintf(fort," %10f",fr->t);
1873 for(i=0; i<norsel; i++)
1874 fprintf(fort," %g",fr->block[enx_i][orsel[i]]);
1875 fprintf(fort,"\n");
1876 }
1877 if (bODT) {
1878 fprintf(fodt," %10f",fr->t);
1879 for(i=0; i<norsel; i++)
1880 fprintf(fodt," %g",fr->block[enx_i][orsel[i]]-oobs[orsel[i]]);
1881 fprintf(fodt,"\n");
1882 }
1883 norfr++;
1884 }
1885 if (bOTEN && fr->nblock>enxORT) {
1886 if (fr->nr[enxORT] != nex*12)
1887 gmx_fatal(FARGS,"Number of orientation experiments in energy file (%g) does not match with the topology (%d)",fr->nr[enxORT]/12,nex);
1888 fprintf(foten," %10f",fr->t);
1889 for(i=0; i<nex; i++)
1890 for(j=0; j<(bOvec?12:3); j++)
1891 fprintf(foten," %g",fr->block[enxORT][i*12+j]);
1892 fprintf(foten,"\n");
1893 }
1894 }
1895 }
1896 }
1897 } while (bCont && (timecheck == 0));
1898
1899 fprintf(stderr,"\n");
1900 close_enx(fp);
1901
1902 ffclose(out);
1903
1904 if (bDRAll)
1905 ffclose(fp_pairs);
1906
1907 if (bORT)
1908 ffclose(fort);
1909 if (bODT)
1910 ffclose(fodt);
1911 if (bORA) {
1912 out = xvgropen(opt2fn("-ora",NFILE,fnm),
1913 "Average calculated orientations",
1914 "Restraint label","",oenv);
1915 if (bOrinst)
1916 fprintf(out,"%s",orinst_sub);
1917 for(i=0; i<nor; i++)
1918 fprintf(out,"%5d %g\n",or_label[i],orient[i]/norfr);
1919 ffclose(out);
1920 }
1921 if (bODA) {
1922 out = xvgropen(opt2fn("-oda",NFILE,fnm),
1923 "Average restraint deviation",
1924 "Restraint label","",oenv);
1925 if (bOrinst)
1926 fprintf(out,"%s",orinst_sub);
1927 for(i=0; i<nor; i++)
1928 fprintf(out,"%5d %g\n",or_label[i],orient[i]/norfr-oobs[i]);
1929 ffclose(out);
1930 }
1931 if (bODR) {
1932 out = xvgropen(opt2fn("-odr",NFILE,fnm),
1933 "RMS orientation restraint deviations",
1934 "Restraint label","",oenv);
1935 if (bOrinst)
1936 fprintf(out,"%s",orinst_sub);
1937 for(i=0; i<nor; i++)
1938 fprintf(out,"%5d %g\n",or_label[i],sqrt(odrms[i]/norfr));
1939 ffclose(out);
1940 }
1941 if (bOTEN)
1942 ffclose(foten);
1943
1944 if (bDisRe) {
1945 analyse_disre(opt2fn("-viol",NFILE,fnm),
1946 teller_disre,violaver,bounds,index,pair,nbounds,oenv);
1947 } else {
1948 analyse_ener(opt2bSet("-corr",NFILE,fnm),opt2fn("-corr",NFILE,fnm),
1949 bFee,bSum,bFluct,opt2parg_bSet("-nmol",npargs,ppa),
1950 bVisco,opt2fn("-vis",NFILE,fnm),
1951 nmol,nconstr,start_step,start_t,frame[cur].step,frame[cur].t,
1952 time,reftemp,&edat,
1953 nset,set,bIsEner,leg,enm,Vaver,ezero,nbmin,nbmax,
1954 oenv);
1955 }
1956 if (opt2bSet("-f2",NFILE,fnm)) {
1957 fec(opt2fn("-f2",NFILE,fnm), opt2fn("-ravg",NFILE,fnm),
1958 reftemp, nset, set, leg, &edat, time ,oenv);
1959 }
1960
1961 {
1962 const char *nxy = "-nxy";
1963
1964 do_view(oenv,opt2fn("-o",NFILE,fnm),nxy);
1965 do_view(oenv,opt2fn_null("-ravg",NFILE,fnm),nxy);
1966 do_view(oenv,opt2fn_null("-ora",NFILE,fnm),nxy);
1967 do_view(oenv,opt2fn_null("-ort",NFILE,fnm),nxy);
1968 do_view(oenv,opt2fn_null("-oda",NFILE,fnm),nxy);
1969 do_view(oenv,opt2fn_null("-odr",NFILE,fnm),nxy);
1970 do_view(oenv,opt2fn_null("-odt",NFILE,fnm),nxy);
1971 do_view(oenv,opt2fn_null("-oten",NFILE,fnm),nxy);
1972 }
1973 thanx(stderr);
1974
1975 return 0;
1976 }
Atempt to write QMMM interface for Gamess-US