?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Convert gmx_mtop_t to C++
[gromacs.git] / src / gromacs / gmxana / gmx_energy.cpp
blob77de9a322276b520d18fd9f480a8495f9fab2745
1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
10 *
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
15 *
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
20 *
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 *
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
33 *
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
36 */
37 #include "gmxpre.h"
38
39 #include <cmath>
40 #include <cstdlib>
41 #include <cstring>
42
43 #include <algorithm>
44
45 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/commandline/viewit.h"
47 #include "gromacs/correlationfunctions/autocorr.h"
48 #include "gromacs/fileio/enxio.h"
49 #include "gromacs/fileio/gmxfio.h"
50 #include "gromacs/fileio/tpxio.h"
51 #include "gromacs/fileio/trxio.h"
52 #include "gromacs/fileio/xvgr.h"
53 #include "gromacs/gmxana/gmx_ana.h"
54 #include "gromacs/gmxana/gstat.h"
55 #include "gromacs/math/functions.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/mdlib/mdebin.h"
59 #include "gromacs/mdtypes/inputrec.h"
60 #include "gromacs/mdtypes/md_enums.h"
61 #include "gromacs/topology/ifunc.h"
62 #include "gromacs/topology/mtop_lookup.h"
63 #include "gromacs/topology/mtop_util.h"
64 #include "gromacs/topology/topology.h"
65 #include "gromacs/utility/arraysize.h"
66 #include "gromacs/utility/cstringutil.h"
67 #include "gromacs/utility/fatalerror.h"
68 #include "gromacs/utility/gmxassert.h"
69 #include "gromacs/utility/pleasecite.h"
70 #include "gromacs/utility/smalloc.h"
71 #include "gromacs/utility/strconvert.h"
72
73 static const int NOTSET = -23451;
74
75 typedef struct {
76 real sum;
77 real sum2;
78 } exactsum_t;
79
80 typedef struct {
81 real *ener;
82 exactsum_t *es;
83 gmx_bool bExactStat;
84 double av;
85 double rmsd;
86 double ee;
87 double slope;
88 } enerdat_t;
89
90 typedef struct {
91 gmx_int64_t nsteps;
92 gmx_int64_t npoints;
93 int nframes;
94 int *step;
95 int *steps;
96 int *points;
97 enerdat_t *s;
98 gmx_bool bHaveSums;
99 } enerdata_t;
100
101 static void done_enerdata_t(int nset, enerdata_t *edat)
102 {
103 sfree(edat->step);
104 sfree(edat->steps);
105 sfree(edat->points);
106 for (int i = 0; i < nset; i++)
107 {
108 sfree(edat->s[i].ener);
109 sfree(edat->s[i].es);
110 }
111 sfree(edat->s);
112 }
113
114 static void chomp(char *buf)
115 {
116 int len = std::strlen(buf);
117
118 while ((len > 0) && (buf[len-1] == '\n'))
119 {
120 buf[len-1] = '\0';
121 len--;
122 }
123 }
124
125 static int *select_by_name(int nre, gmx_enxnm_t *nm, int *nset)
126 {
127 gmx_bool *bE;
128 int k, kk, j, i, nmatch, nind, nss;
129 int *set;
130 gmx_bool bEOF, bVerbose = TRUE, bLong = FALSE;
131 char *ptr, buf[STRLEN];
132 const char *fm4 = "%3d %-14s";
133 const char *fm2 = "%3d %-34s";
134 char **newnm = nullptr;
135
136 if ((getenv("GMX_ENER_VERBOSE")) != nullptr)
137 {
138 bVerbose = FALSE;
139 }
140
141 fprintf(stderr, "\n");
142 fprintf(stderr, "Select the terms you want from the following list by\n");
143 fprintf(stderr, "selecting either (part of) the name or the number or a combination.\n");
144 fprintf(stderr, "End your selection with an empty line or a zero.\n");
145 fprintf(stderr, "-------------------------------------------------------------------\n");
146
147 snew(newnm, nre);
148 j = 0;
149 for (k = 0; k < nre; k++)
150 {
151 newnm[k] = gmx_strdup(nm[k].name);
152 /* Insert dashes in all the names */
153 while ((ptr = std::strchr(newnm[k], ' ')) != nullptr)
154 {
155 *ptr = '-';
156 }
157 if (bVerbose)
158 {
159 if (j == 0)
160 {
161 if (k > 0)
162 {
163 fprintf(stderr, "\n");
164 }
165 bLong = FALSE;
166 for (kk = k; kk < k+4; kk++)
167 {
168 if (kk < nre && std::strlen(nm[kk].name) > 14)
169 {
170 bLong = TRUE;
171 }
172 }
173 }
174 else
175 {
176 fprintf(stderr, " ");
177 }
178 if (!bLong)
179 {
180 fprintf(stderr, fm4, k+1, newnm[k]);
181 j++;
182 if (j == 4)
183 {
184 j = 0;
185 }
186 }
187 else
188 {
189 fprintf(stderr, fm2, k+1, newnm[k]);
190 j++;
191 if (j == 2)
192 {
193 j = 0;
194 }
195 }
196 }
197 }
198 if (bVerbose)
199 {
200 fprintf(stderr, "\n\n");
201 }
202
203 snew(bE, nre);
204
205 bEOF = FALSE;
206 while (!bEOF && (fgets2(buf, STRLEN-1, stdin)))
207 {
208 /* Remove newlines */
209 chomp(buf);
210
211 /* Remove spaces */
212 trim(buf);
213
214 /* Empty line means end of input */
215 bEOF = (std::strlen(buf) == 0);
216 if (!bEOF)
217 {
218 ptr = buf;
219 do
220 {
221 if (!bEOF)
222 {
223 /* First try to read an integer */
224 nss = sscanf(ptr, "%d", &nind);
225 if (nss == 1)
226 {
227 /* Zero means end of input */
228 if (nind == 0)
229 {
230 bEOF = TRUE;
231 }
232 else if ((1 <= nind) && (nind <= nre))
233 {
234 bE[nind-1] = TRUE;
235 }
236 else
237 {
238 fprintf(stderr, "number %d is out of range\n", nind);
239 }
240 }
241 else
242 {
243 /* Now try to read a string */
244 nmatch = 0;
245 for (nind = 0; nind < nre; nind++)
246 {
247 if (gmx_strcasecmp(newnm[nind], ptr) == 0)
248 {
249 bE[nind] = TRUE;
250 nmatch++;
251 }
252 }
253 if (nmatch == 0)
254 {
255 i = std::strlen(ptr);
256 nmatch = 0;
257 for (nind = 0; nind < nre; nind++)
258 {
259 if (gmx_strncasecmp(newnm[nind], ptr, i) == 0)
260 {
261 bE[nind] = TRUE;
262 nmatch++;
263 }
264 }
265 if (nmatch == 0)
266 {
267 fprintf(stderr, "String '%s' does not match anything\n", ptr);
268 }
269 }
270 }
271 }
272 /* Look for the first space, and remove spaces from there */
273 if ((ptr = std::strchr(ptr, ' ')) != nullptr)
274 {
275 trim(ptr);
276 }
277 }
278 while (!bEOF && (ptr && (std::strlen(ptr) > 0)));
279 }
280 }
281
282 snew(set, nre);
283 for (i = (*nset) = 0; (i < nre); i++)
284 {
285 if (bE[i])
286 {
287 set[(*nset)++] = i;
288 }
289 }
290
291 sfree(bE);
292
293 if (*nset == 0)
294 {
295 gmx_fatal(FARGS, "No energy terms selected");
296 }
297
298 for (i = 0; (i < nre); i++)
299 {
300 sfree(newnm[i]);
301 }
302 sfree(newnm);
303
304 return set;
305 }
306
307 static void get_dhdl_parms(const char *topnm, t_inputrec *ir)
308 {
309 gmx_mtop_t mtop;
310 int natoms;
311 matrix box;
312
313 /* all we need is the ir to be able to write the label */
314 read_tpx(topnm, ir, box, &natoms, nullptr, nullptr, &mtop);
315 }
316
317 static void einstein_visco(const char *fn, const char *fni, int nsets,
318 int nint, real **eneint,
319 real V, real T, double dt,
320 const gmx_output_env_t *oenv)
321 {
322 FILE *fp0, *fp1;
323 double av[4], avold[4];
324 double fac, di;
325 int i, j, m, nf4;
326
327 nf4 = nint/4 + 1;
328
329 for (i = 0; i <= nsets; i++)
330 {
331 avold[i] = 0;
332 }
333 fp0 = xvgropen(fni, "Shear viscosity integral",
334 "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N ps)", oenv);
335 fp1 = xvgropen(fn, "Shear viscosity using Einstein relation",
336 "Time (ps)", "(kg m\\S-1\\N s\\S-1\\N)", oenv);
337 for (i = 0; i < nf4; i++)
338 {
339 for (m = 0; m <= nsets; m++)
340 {
341 av[m] = 0;
342 }
343 for (j = 0; j < nint - i; j++)
344 {
345 for (m = 0; m < nsets; m++)
346 {
347 di = gmx::square(eneint[m][j+i] - eneint[m][j]);
348
349 av[m] += di;
350 av[nsets] += di/nsets;
351 }
352 }
353 /* Convert to SI for the viscosity */
354 fac = (V*NANO*NANO*NANO*PICO*1e10)/(2*BOLTZMANN*T)/(nint - i);
355 fprintf(fp0, "%10g", i*dt);
356 for (m = 0; (m <= nsets); m++)
357 {
358 av[m] = fac*av[m];
359 fprintf(fp0, " %10g", av[m]);
360 }
361 fprintf(fp0, "\n");
362 fprintf(fp1, "%10g", (i + 0.5)*dt);
363 for (m = 0; (m <= nsets); m++)
364 {
365 fprintf(fp1, " %10g", (av[m]-avold[m])/dt);
366 avold[m] = av[m];
367 }
368 fprintf(fp1, "\n");
369 }
370 xvgrclose(fp0);
371 xvgrclose(fp1);
372 }
373
374 typedef struct {
375 gmx_int64_t np;
376 double sum;
377 double sav;
378 double sav2;
379 } ee_sum_t;
380
381 typedef struct {
382 int b;
383 ee_sum_t sum;
384 gmx_int64_t nst;
385 gmx_int64_t nst_min;
386 } ener_ee_t;
387
388 static void clear_ee_sum(ee_sum_t *ees)
389 {
390 ees->sav = 0;
391 ees->sav2 = 0;
392 ees->np = 0;
393 ees->sum = 0;
394 }
395
396 static void add_ee_sum(ee_sum_t *ees, double sum, int np)
397 {
398 ees->np += np;
399 ees->sum += sum;
400 }
401
402 static void add_ee_av(ee_sum_t *ees)
403 {
404 double av;
405
406 av = ees->sum/ees->np;
407 ees->sav += av;
408 ees->sav2 += av*av;
409 ees->np = 0;
410 ees->sum = 0;
411 }
412
413 static double calc_ee2(int nb, ee_sum_t *ees)
414 {
415 return (ees->sav2/nb - gmx::square(ees->sav/nb))/(nb - 1);
416 }
417
418 static void set_ee_av(ener_ee_t *eee)
419 {
420 if (debug)
421 {
422 char buf[STEPSTRSIZE];
423 fprintf(debug, "Storing average for err.est.: %s steps\n",
424 gmx_step_str(eee->nst, buf));
425 }
426 add_ee_av(&eee->sum);
427 eee->b++;
428 if (eee->b == 1 || eee->nst < eee->nst_min)
429 {
430 eee->nst_min = eee->nst;
431 }
432 eee->nst = 0;
433 }
434
435 static void calc_averages(int nset, enerdata_t *edat, int nbmin, int nbmax)
436 {
437 int nb, i, f, nee;
438 double sum, sum2, sump, see2;
439 gmx_int64_t np, p, bound_nb;
440 enerdat_t *ed;
441 exactsum_t *es;
442 gmx_bool bAllZero;
443 double x, sx, sy, sxx, sxy;
444 ener_ee_t *eee;
445
446 /* Check if we have exact statistics over all points */
447 for (i = 0; i < nset; i++)
448 {
449 ed = &edat->s[i];
450 ed->bExactStat = FALSE;
451 if (edat->bHaveSums)
452 {
453 /* All energy file sum entries 0 signals no exact sums.
454 * But if all energy values are 0, we still have exact sums.
455 */
456 bAllZero = TRUE;
457 for (f = 0; f < edat->nframes && !ed->bExactStat; f++)
458 {
459 if (ed->ener[i] != 0)
460 {
461 bAllZero = FALSE;
462 }
463 ed->bExactStat = (ed->es[f].sum != 0);
464 }
465 if (bAllZero)
466 {
467 ed->bExactStat = TRUE;
468 }
469 }
470 }
471
472 snew(eee, nbmax+1);
473 for (i = 0; i < nset; i++)
474 {
475 ed = &edat->s[i];
476
477 sum = 0;
478 sum2 = 0;
479 np = 0;
480 sx = 0;
481 sy = 0;
482 sxx = 0;
483 sxy = 0;
484 for (nb = nbmin; nb <= nbmax; nb++)
485 {
486 eee[nb].b = 0;
487 clear_ee_sum(&eee[nb].sum);
488 eee[nb].nst = 0;
489 eee[nb].nst_min = 0;
490 }
491 for (f = 0; f < edat->nframes; f++)
492 {
493 es = &ed->es[f];
494
495 if (ed->bExactStat)
496 {
497 /* Add the sum and the sum of variances to the totals. */
498 p = edat->points[f];
499 sump = es->sum;
500 sum2 += es->sum2;
501 if (np > 0)
502 {
503 sum2 += gmx::square(sum/np - (sum + es->sum)/(np + p))
504 *np*(np + p)/p;
505 }
506 }
507 else
508 {
509 /* Add a single value to the sum and sum of squares. */
510 p = 1;
511 sump = ed->ener[f];
512 sum2 += gmx::square(sump);
513 }
514
515 /* sum has to be increased after sum2 */
516 np += p;
517 sum += sump;
518
519 /* For the linear regression use variance 1/p.
520 * Note that sump is the sum, not the average, so we don't need p*.
521 */
522 x = edat->step[f] - 0.5*(edat->steps[f] - 1);
523 sx += p*x;
524 sy += sump;
525 sxx += p*x*x;
526 sxy += x*sump;
527
528 for (nb = nbmin; nb <= nbmax; nb++)
529 {
530 /* Check if the current end step is closer to the desired
531 * block boundary than the next end step.
532 */
533 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
534 if (eee[nb].nst > 0 &&
535 bound_nb - edat->step[f-1]*nb < edat->step[f]*nb - bound_nb)
536 {
537 set_ee_av(&eee[nb]);
538 }
539 if (f == 0)
540 {
541 eee[nb].nst = 1;
542 }
543 else
544 {
545 eee[nb].nst += edat->step[f] - edat->step[f-1];
546 }
547 if (ed->bExactStat)
548 {
549 add_ee_sum(&eee[nb].sum, es->sum, edat->points[f]);
550 }
551 else
552 {
553 add_ee_sum(&eee[nb].sum, edat->s[i].ener[f], 1);
554 }
555 bound_nb = (edat->step[0]-1)*nb + edat->nsteps*(eee[nb].b+1);
556 if (edat->step[f]*nb >= bound_nb)
557 {
558 set_ee_av(&eee[nb]);
559 }
560 }
561 }
562
563 edat->s[i].av = sum/np;
564 if (ed->bExactStat)
565 {
566 edat->s[i].rmsd = std::sqrt(sum2/np);
567 }
568 else
569 {
570 edat->s[i].rmsd = std::sqrt(sum2/np - gmx::square(edat->s[i].av));
571 }
572
573 if (edat->nframes > 1)
574 {
575 edat->s[i].slope = (np*sxy - sx*sy)/(np*sxx - sx*sx);
576 }
577 else
578 {
579 edat->s[i].slope = 0;
580 }
581
582 nee = 0;
583 see2 = 0;
584 for (nb = nbmin; nb <= nbmax; nb++)
585 {
586 /* Check if we actually got nb blocks and if the smallest
587 * block is not shorter than 80% of the average.
588 */
589 if (debug)
590 {
591 char buf1[STEPSTRSIZE], buf2[STEPSTRSIZE];
592 fprintf(debug, "Requested %d blocks, we have %d blocks, min %s nsteps %s\n",
593 nb, eee[nb].b,
594 gmx_step_str(eee[nb].nst_min, buf1),
595 gmx_step_str(edat->nsteps, buf2));
596 }
597 if (eee[nb].b == nb && 5*nb*eee[nb].nst_min >= 4*edat->nsteps)
598 {
599 see2 += calc_ee2(nb, &eee[nb].sum);
600 nee++;
601 }
602 }
603 if (nee > 0)
604 {
605 edat->s[i].ee = std::sqrt(see2/nee);
606 }
607 else
608 {
609 edat->s[i].ee = -1;
610 }
611 }
612 sfree(eee);
613 }
614
615 static enerdata_t *calc_sum(int nset, enerdata_t *edat, int nbmin, int nbmax)
616 {
617 enerdata_t *esum;
618 enerdat_t *s;
619 int f, i;
620 double sum;
621
622 snew(esum, 1);
623 *esum = *edat;
624 snew(esum->s, 1);
625 s = &esum->s[0];
626 snew(s->ener, esum->nframes);
627 snew(s->es, esum->nframes);
628
629 s->bExactStat = TRUE;
630 s->slope = 0;
631 for (i = 0; i < nset; i++)
632 {
633 if (!edat->s[i].bExactStat)
634 {
635 s->bExactStat = FALSE;
636 }
637 s->slope += edat->s[i].slope;
638 }
639
640 for (f = 0; f < edat->nframes; f++)
641 {
642 sum = 0;
643 for (i = 0; i < nset; i++)
644 {
645 sum += edat->s[i].ener[f];
646 }
647 s->ener[f] = sum;
648 sum = 0;
649 for (i = 0; i < nset; i++)
650 {
651 sum += edat->s[i].es[f].sum;
652 }
653 s->es[f].sum = sum;
654 s->es[f].sum2 = 0;
655 }
656
657 calc_averages(1, esum, nbmin, nbmax);
658
659 return esum;
660 }
661
662 static void ee_pr(double ee, int buflen, char *buf)
663 {
664 snprintf(buf, buflen, "%s", "--");
665 if (ee >= 0)
666 {
667 /* Round to two decimals by printing. */
668 char tmp[100];
669 snprintf(tmp, sizeof(tmp), "%.1e", ee);
670 double rnd = gmx::doubleFromString(tmp);
671 snprintf(buf, buflen, "%g", rnd);
672 }
673 }
674
675 static void remove_drift(int nset, int nbmin, int nbmax, real dt, enerdata_t *edat)
676 {
677 /* Remove the drift by performing a fit to y = ax+b.
678 Uses 5 iterations. */
679 int i, j, k;
680 double delta;
681
682 edat->npoints = edat->nframes;
683 edat->nsteps = edat->nframes;
684
685 for (k = 0; (k < 5); k++)
686 {
687 for (i = 0; (i < nset); i++)
688 {
689 delta = edat->s[i].slope*dt;
690
691 if (nullptr != debug)
692 {
693 fprintf(debug, "slope for set %d is %g\n", i, edat->s[i].slope);
694 }
695
696 for (j = 0; (j < edat->nframes); j++)
697 {
698 edat->s[i].ener[j] -= j*delta;
699 edat->s[i].es[j].sum = 0;
700 edat->s[i].es[j].sum2 = 0;
701 }
702 }
703 calc_averages(nset, edat, nbmin, nbmax);
704 }
705 }
706
707 static void calc_fluctuation_props(FILE *fp,
708 gmx_bool bDriftCorr, real dt,
709 int nset, int nmol,
710 char **leg, enerdata_t *edat,
711 int nbmin, int nbmax)
712 {
713 int i, j;
714 double vv, v, h, varv, hh, varh, tt, cv, cp, alpha, kappa, dcp, varet;
715 double NANO3;
716 enum {
717 eVol, eEnth, eTemp, eEtot, eNR
718 };
719 const char *my_ener[] = { "Volume", "Enthalpy", "Temperature", "Total Energy" };
720 int ii[eNR];
721
722 NANO3 = NANO*NANO*NANO;
723 if (!bDriftCorr)
724 {
725 fprintf(fp, "\nYou may want to use the -driftcorr flag in order to correct\n"
726 "for spurious drift in the graphs. Note that this is not\n"
727 "a substitute for proper equilibration and sampling!\n");
728 }
729 else
730 {
731 remove_drift(nset, nbmin, nbmax, dt, edat);
732 }
733 for (i = 0; (i < eNR); i++)
734 {
735 for (ii[i] = 0; (ii[i] < nset &&
736 (gmx_strcasecmp(leg[ii[i]], my_ener[i]) != 0)); ii[i]++)
737 {
738 ;
739 }
740 /* if (ii[i] < nset)
741 fprintf(fp,"Found %s data.\n",my_ener[i]);
742 */ }
743 /* Compute it all! */
744 alpha = kappa = cp = dcp = cv = NOTSET;
745
746 /* Temperature */
747 tt = NOTSET;
748 if (ii[eTemp] < nset)
749 {
750 tt = edat->s[ii[eTemp]].av;
751 }
752 /* Volume */
753 vv = varv = NOTSET;
754 if ((ii[eVol] < nset) && (ii[eTemp] < nset))
755 {
756 vv = edat->s[ii[eVol]].av*NANO3;
757 varv = gmx::square(edat->s[ii[eVol]].rmsd*NANO3);
758 kappa = (varv/vv)/(BOLTZMANN*tt);
759 }
760 /* Enthalpy */
761 hh = varh = NOTSET;
762 if ((ii[eEnth] < nset) && (ii[eTemp] < nset))
763 {
764 hh = KILO*edat->s[ii[eEnth]].av/AVOGADRO;
765 varh = gmx::square(KILO*edat->s[ii[eEnth]].rmsd/AVOGADRO);
766 cp = AVOGADRO*((varh/nmol)/(BOLTZMANN*tt*tt));
767 }
768 /* Total energy */
769 if ((ii[eEtot] < nset) && (hh == NOTSET) && (tt != NOTSET))
770 {
771 /* Only compute cv in constant volume runs, which we can test
772 by checking whether the enthalpy was computed.
773 */
774 varet = gmx::square(edat->s[ii[eEtot]].rmsd);
775 cv = KILO*((varet/nmol)/(BOLTZ*tt*tt));
776 }
777 /* Alpha, dcp */
778 if ((ii[eVol] < nset) && (ii[eEnth] < nset) && (ii[eTemp] < nset))
779 {
780 double v_sum, h_sum, vh_sum, v_aver, h_aver, vh_aver;
781 vh_sum = v_sum = h_sum = 0;
782 for (j = 0; (j < edat->nframes); j++)
783 {
784 v = edat->s[ii[eVol]].ener[j]*NANO3;
785 h = KILO*edat->s[ii[eEnth]].ener[j]/AVOGADRO;
786 v_sum += v;
787 h_sum += h;
788 vh_sum += (v*h);
789 }
790 vh_aver = vh_sum / edat->nframes;
791 v_aver = v_sum / edat->nframes;
792 h_aver = h_sum / edat->nframes;
793 alpha = (vh_aver-v_aver*h_aver)/(v_aver*BOLTZMANN*tt*tt);
794 dcp = (v_aver*AVOGADRO/nmol)*tt*gmx::square(alpha)/(kappa);
795 }
796
797 if (tt != NOTSET)
798 {
799 if (nmol < 2)
800 {
801 fprintf(fp, "\nWARNING: nmol = %d, this may not be what you want.\n",
802 nmol);
803 }
804 fprintf(fp, "\nTemperature dependent fluctuation properties at T = %g.\n", tt);
805 fprintf(fp, "\nHeat capacities obtained from fluctuations do *not* include\n");
806 fprintf(fp, "quantum corrections. If you want to get a more accurate estimate\n");
807 fprintf(fp, "please use the g_dos program.\n\n");
808 fprintf(fp, "WARNING: Please verify that your simulations are converged and perform\n"
809 "a block-averaging error analysis (not implemented in g_energy yet)\n");
810
811 if (debug != nullptr)
812 {
813 if (varv != NOTSET)
814 {
815 fprintf(fp, "varv = %10g (m^6)\n", varv*AVOGADRO/nmol);
816 }
817 }
818 if (vv != NOTSET)
819 {
820 fprintf(fp, "Volume = %10g m^3/mol\n",
821 vv*AVOGADRO/nmol);
822 }
823 if (varh != NOTSET)
824 {
825 fprintf(fp, "Enthalpy = %10g kJ/mol\n",
826 hh*AVOGADRO/(KILO*nmol));
827 }
828 if (alpha != NOTSET)
829 {
830 fprintf(fp, "Coefficient of Thermal Expansion Alpha_P = %10g (1/K)\n",
831 alpha);
832 }
833 if (kappa != NOTSET)
834 {
835 fprintf(fp, "Isothermal Compressibility Kappa = %10g (m^3/J)\n",
836 kappa);
837 fprintf(fp, "Adiabatic bulk modulus = %10g (J/m^3)\n",
838 1.0/kappa);
839 }
840 if (cp != NOTSET)
841 {
842 fprintf(fp, "Heat capacity at constant pressure Cp = %10g J/(mol K)\n",
843 cp);
844 }
845 if (cv != NOTSET)
846 {
847 fprintf(fp, "Heat capacity at constant volume Cv = %10g J/(mol K)\n",
848 cv);
849 }
850 if (dcp != NOTSET)
851 {
852 fprintf(fp, "Cp-Cv = %10g J/(mol K)\n",
853 dcp);
854 }
855 please_cite(fp, "Allen1987a");
856 }
857 else
858 {
859 fprintf(fp, "You should select the temperature in order to obtain fluctuation properties.\n");
860 }
861 }
862
863 static void analyse_ener(gmx_bool bCorr, const char *corrfn,
864 const char *eviscofn, const char *eviscoifn,
865 gmx_bool bFee, gmx_bool bSum, gmx_bool bFluct,
866 gmx_bool bVisco, const char *visfn, int nmol,
867 gmx_int64_t start_step, double start_t,
868 gmx_int64_t step, double t,
869 real reftemp,
870 enerdata_t *edat,
871 int nset, int set[], gmx_bool *bIsEner,
872 char **leg, gmx_enxnm_t *enm,
873 real Vaver, real ezero,
874 int nbmin, int nbmax,
875 const gmx_output_env_t *oenv)
876 {
877 FILE *fp;
878 /* Check out the printed manual for equations! */
879 double Dt, aver, stddev, errest, delta_t, totaldrift;
880 enerdata_t *esum = nullptr;
881 real integral, intBulk, Temp = 0, Pres = 0;
882 real pr_aver, pr_stddev, pr_errest;
883 double beta = 0, expE, expEtot, *fee = nullptr;
884 gmx_int64_t nsteps;
885 int nexact, nnotexact;
886 int i, j, nout;
887 char buf[256], eebuf[100];
888
889 nsteps = step - start_step + 1;
890 if (nsteps < 1)
891 {
892 fprintf(stdout, "Not enough steps (%s) for statistics\n",
893 gmx_step_str(nsteps, buf));
894 }
895 else
896 {
897 /* Calculate the time difference */
898 delta_t = t - start_t;
899
900 fprintf(stdout, "\nStatistics over %s steps [ %.4f through %.4f ps ], %d data sets\n",
901 gmx_step_str(nsteps, buf), start_t, t, nset);
902
903 calc_averages(nset, edat, nbmin, nbmax);
904
905 if (bSum)
906 {
907 esum = calc_sum(nset, edat, nbmin, nbmax);
908 }
909
910 if (!edat->bHaveSums)
911 {
912 nexact = 0;
913 nnotexact = nset;
914 }
915 else
916 {
917 nexact = 0;
918 nnotexact = 0;
919 for (i = 0; (i < nset); i++)
920 {
921 if (edat->s[i].bExactStat)
922 {
923 nexact++;
924 }
925 else
926 {
927 nnotexact++;
928 }
929 }
930 }
931
932 if (nnotexact == 0)
933 {
934 fprintf(stdout, "All statistics are over %s points\n",
935 gmx_step_str(edat->npoints, buf));
936 }
937 else if (nexact == 0 || edat->npoints == edat->nframes)
938 {
939 fprintf(stdout, "All statistics are over %d points (frames)\n",
940 edat->nframes);
941 }
942 else
943 {
944 fprintf(stdout, "The term%s", nnotexact == 1 ? "" : "s");
945 for (i = 0; (i < nset); i++)
946 {
947 if (!edat->s[i].bExactStat)
948 {
949 fprintf(stdout, " '%s'", leg[i]);
950 }
951 }
952 fprintf(stdout, " %s has statistics over %d points (frames)\n",
953 nnotexact == 1 ? "is" : "are", edat->nframes);
954 fprintf(stdout, "All other statistics are over %s points\n",
955 gmx_step_str(edat->npoints, buf));
956 }
957 fprintf(stdout, "\n");
958
959 fprintf(stdout, "%-24s %10s %10s %10s %10s",
960 "Energy", "Average", "Err.Est.", "RMSD", "Tot-Drift");
961 if (bFee)
962 {
963 fprintf(stdout, " %10s\n", "-kT ln<e^(E/kT)>");
964 }
965 else
966 {
967 fprintf(stdout, "\n");
968 }
969 fprintf(stdout, "-------------------------------------------------------------------------------\n");
970
971 /* Initiate locals, only used with -sum */
972 expEtot = 0;
973 if (bFee)
974 {
975 beta = 1.0/(BOLTZ*reftemp);
976 snew(fee, nset);
977 }
978 for (i = 0; (i < nset); i++)
979 {
980 aver = edat->s[i].av;
981 stddev = edat->s[i].rmsd;
982 errest = edat->s[i].ee;
983
984 if (bFee)
985 {
986 expE = 0;
987 for (j = 0; (j < edat->nframes); j++)
988 {
989 expE += std::exp(beta*(edat->s[i].ener[j] - aver)/nmol);
990 }
991 if (bSum)
992 {
993 expEtot += expE/edat->nframes;
994 }
995
996 fee[i] = std::log(expE/edat->nframes)/beta + aver/nmol;
997 }
998 if (std::strstr(leg[i], "empera") != nullptr)
999 {
1000 Temp = aver;
1001 }
1002 else if (std::strstr(leg[i], "olum") != nullptr)
1003 {
1004 Vaver = aver;
1005 }
1006 else if (std::strstr(leg[i], "essure") != nullptr)
1007 {
1008 Pres = aver;
1009 }
1010 if (bIsEner[i])
1011 {
1012 pr_aver = aver/nmol-ezero;
1013 pr_stddev = stddev/nmol;
1014 pr_errest = errest/nmol;
1015 }
1016 else
1017 {
1018 pr_aver = aver;
1019 pr_stddev = stddev;
1020 pr_errest = errest;
1021 }
1022
1023 /* Multiply the slope in steps with the number of steps taken */
1024 totaldrift = (edat->nsteps - 1)*edat->s[i].slope;
1025 if (bIsEner[i])
1026 {
1027 totaldrift /= nmol;
1028 }
1029
1030 ee_pr(pr_errest, sizeof(eebuf), eebuf);
1031 fprintf(stdout, "%-24s %10g %10s %10g %10g",
1032 leg[i], pr_aver, eebuf, pr_stddev, totaldrift);
1033 if (bFee)
1034 {
1035 fprintf(stdout, " %10g", fee[i]);
1036 }
1037
1038 fprintf(stdout, " (%s)\n", enm[set[i]].unit);
1039
1040 if (bFluct)
1041 {
1042 for (j = 0; (j < edat->nframes); j++)
1043 {
1044 edat->s[i].ener[j] -= aver;
1045 }
1046 }
1047 }
1048 if (bSum)
1049 {
1050 totaldrift = (edat->nsteps - 1)*esum->s[0].slope;
1051 ee_pr(esum->s[0].ee/nmol, sizeof(eebuf), eebuf);
1052 fprintf(stdout, "%-24s %10g %10s %10s %10g (%s)",
1053 "Total", esum->s[0].av/nmol, eebuf,
1054 "--", totaldrift/nmol, enm[set[0]].unit);
1055 /* pr_aver,pr_stddev,a,totaldrift */
1056 if (bFee)
1057 {
1058 fprintf(stdout, " %10g %10g\n",
1059 std::log(expEtot)/beta + esum->s[0].av/nmol, std::log(expEtot)/beta);
1060 }
1061 else
1062 {
1063 fprintf(stdout, "\n");
1064 }
1065 }
1066
1067 /* Do correlation function */
1068 if (edat->nframes > 1)
1069 {
1070 Dt = delta_t/(edat->nframes - 1);
1071 }
1072 else
1073 {
1074 Dt = 0;
1075 }
1076 if (bVisco)
1077 {
1078 const char* leg[] = { "Shear", "Bulk" };
1079 real factor;
1080 real **eneset;
1081 real **eneint;
1082
1083 /* Assume pressure tensor is in Pxx Pxy Pxz Pyx Pyy Pyz Pzx Pzy Pzz */
1084
1085 /* Symmetrise tensor! (and store in first three elements)
1086 * And subtract average pressure!
1087 */
1088 snew(eneset, 12);
1089 for (i = 0; i < 12; i++)
1090 {
1091 snew(eneset[i], edat->nframes);
1092 }
1093 for (i = 0; (i < edat->nframes); i++)
1094 {
1095 eneset[0][i] = 0.5*(edat->s[1].ener[i]+edat->s[3].ener[i]);
1096 eneset[1][i] = 0.5*(edat->s[2].ener[i]+edat->s[6].ener[i]);
1097 eneset[2][i] = 0.5*(edat->s[5].ener[i]+edat->s[7].ener[i]);
1098 for (j = 3; j <= 11; j++)
1099 {
1100 eneset[j][i] = edat->s[j].ener[i];
1101 }
1102 eneset[11][i] -= Pres;
1103 }
1104
1105 /* Determine integrals of the off-diagonal pressure elements */
1106 snew(eneint, 3);
1107 for (i = 0; i < 3; i++)
1108 {
1109 snew(eneint[i], edat->nframes + 1);
1110 }
1111 eneint[0][0] = 0;
1112 eneint[1][0] = 0;
1113 eneint[2][0] = 0;
1114 for (i = 0; i < edat->nframes; i++)
1115 {
1116 eneint[0][i+1] = eneint[0][i] + 0.5*(edat->s[1].es[i].sum + edat->s[3].es[i].sum)*Dt/edat->points[i];
1117 eneint[1][i+1] = eneint[1][i] + 0.5*(edat->s[2].es[i].sum + edat->s[6].es[i].sum)*Dt/edat->points[i];
1118 eneint[2][i+1] = eneint[2][i] + 0.5*(edat->s[5].es[i].sum + edat->s[7].es[i].sum)*Dt/edat->points[i];
1119 }
1120
1121 einstein_visco(eviscofn, eviscoifn,
1122 3, edat->nframes+1, eneint, Vaver, Temp, Dt, oenv);
1123
1124 for (i = 0; i < 3; i++)
1125 {
1126 sfree(eneint[i]);
1127 }
1128 sfree(eneint);
1129
1130 /*do_autocorr(corrfn,buf,nenergy,3,eneset,Dt,eacNormal,TRUE);*/
1131 /* Do it for shear viscosity */
1132 std::strcpy(buf, "Shear Viscosity");
1133 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 3,
1134 (edat->nframes+1)/2, eneset, Dt,
1135 eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1136
1137 /* Now for bulk viscosity */
1138 std::strcpy(buf, "Bulk Viscosity");
1139 low_do_autocorr(corrfn, oenv, buf, edat->nframes, 1,
1140 (edat->nframes+1)/2, &(eneset[11]), Dt,
1141 eacNormal, 1, TRUE, FALSE, FALSE, 0.0, 0.0, 0);
1142
1143 factor = (Vaver*1e-26/(BOLTZMANN*Temp))*Dt;
1144 fp = xvgropen(visfn, buf, "Time (ps)", "\\8h\\4 (cp)", oenv);
1145 xvgr_legend(fp, asize(leg), leg, oenv);
1146
1147 /* Use trapezium rule for integration */
1148 integral = 0;
1149 intBulk = 0;
1150 nout = get_acfnout();
1151 if ((nout < 2) || (nout >= edat->nframes/2))
1152 {
1153 nout = edat->nframes/2;
1154 }
1155 for (i = 1; (i < nout); i++)
1156 {
1157 integral += 0.5*(eneset[0][i-1] + eneset[0][i])*factor;
1158 intBulk += 0.5*(eneset[11][i-1] + eneset[11][i])*factor;
1159 fprintf(fp, "%10g %10g %10g\n", (i*Dt), integral, intBulk);
1160 }
1161 xvgrclose(fp);
1162
1163 for (i = 0; i < 12; i++)
1164 {
1165 sfree(eneset[i]);
1166 }
1167 sfree(eneset);
1168 }
1169 else if (bCorr)
1170 {
1171 if (bFluct)
1172 {
1173 std::strcpy(buf, "Autocorrelation of Energy Fluctuations");
1174 }
1175 else
1176 {
1177 std::strcpy(buf, "Energy Autocorrelation");
1178 }
1179 #if 0
1180 do_autocorr(corrfn, oenv, buf, edat->nframes,
1181 bSum ? 1 : nset,
1182 bSum ? &edat->s[nset-1].ener : eneset,
1183 (delta_t/edat->nframes), eacNormal, FALSE);
1184 #endif
1185 }
1186 }
1187 }
1188
1189 static void print_time(FILE *fp, double t)
1190 {
1191 fprintf(fp, "%12.6f", t);
1192 }
1193
1194 static void print1(FILE *fp, gmx_bool bDp, real e)
1195 {
1196 if (bDp)
1197 {
1198 fprintf(fp, " %16.12f", e);
1199 }
1200 else
1201 {
1202 fprintf(fp, " %10.6f", e);
1203 }
1204 }
1205
1206 static void fec(const char *ene2fn, const char *runavgfn,
1207 real reftemp, int nset, int set[], char *leg[],
1208 enerdata_t *edat, double time[],
1209 const gmx_output_env_t *oenv)
1210 {
1211 const char * ravgleg[] = {
1212 "\\8D\\4E = E\\sB\\N-E\\sA\\N",
1213 "<e\\S-\\8D\\4E/kT\\N>\\s0..t\\N"
1214 };
1215 FILE *fp;
1216 ener_file_t enx;
1217 int timecheck, nenergy, nenergy2, maxenergy;
1218 int i, j;
1219 gmx_bool bCont;
1220 real aver, beta;
1221 real **eneset2;
1222 double dE, sum;
1223 gmx_enxnm_t *enm = nullptr;
1224 t_enxframe *fr;
1225 char buf[22];
1226
1227 /* read second energy file */
1228 snew(fr, 1);
1229 enm = nullptr;
1230 enx = open_enx(ene2fn, "r");
1231 do_enxnms(enx, &(fr->nre), &enm);
1232
1233 snew(eneset2, nset+1);
1234 nenergy2 = 0;
1235 maxenergy = 0;
1236 timecheck = 0;
1237 do
1238 {
1239 /* This loop searches for the first frame (when -b option is given),
1240 * or when this has been found it reads just one energy frame
1241 */
1242 do
1243 {
1244 bCont = do_enx(enx, fr);
1245
1246 if (bCont)
1247 {
1248 timecheck = check_times(fr->t);
1249 }
1250
1251 }
1252 while (bCont && (timecheck < 0));
1253
1254 /* Store energies for analysis afterwards... */
1255 if ((timecheck == 0) && bCont)
1256 {
1257 if (fr->nre > 0)
1258 {
1259 if (nenergy2 >= maxenergy)
1260 {
1261 maxenergy += 1000;
1262 for (i = 0; i <= nset; i++)
1263 {
1264 srenew(eneset2[i], maxenergy);
1265 }
1266 }
1267 GMX_RELEASE_ASSERT(time != nullptr, "trying to dereference NULL time pointer");
1268
1269 if (fr->t != time[nenergy2])
1270 {
1271 fprintf(stderr, "\nWARNING time mismatch %g!=%g at frame %s\n",
1272 fr->t, time[nenergy2], gmx_step_str(fr->step, buf));
1273 }
1274 for (i = 0; i < nset; i++)
1275 {
1276 eneset2[i][nenergy2] = fr->ener[set[i]].e;
1277 }
1278 nenergy2++;
1279 }
1280 }
1281 }
1282 while (bCont && (timecheck == 0));
1283
1284 /* check */
1285 if (edat->nframes != nenergy2)
1286 {
1287 fprintf(stderr, "\nWARNING file length mismatch %d!=%d\n",
1288 edat->nframes, nenergy2);
1289 }
1290 nenergy = std::min(edat->nframes, nenergy2);
1291
1292 /* calculate fe difference dF = -kT ln < exp(-(E_B-E_A)/kT) >_A */
1293 fp = nullptr;
1294 if (runavgfn)
1295 {
1296 fp = xvgropen(runavgfn, "Running average free energy difference",
1297 "Time (" unit_time ")", "\\8D\\4E (" unit_energy ")", oenv);
1298 xvgr_legend(fp, asize(ravgleg), ravgleg, oenv);
1299 }
1300 fprintf(stdout, "\n%-24s %10s\n",
1301 "Energy", "dF = -kT ln < exp(-(EB-EA)/kT) >A");
1302 sum = 0;
1303 beta = 1.0/(BOLTZ*reftemp);
1304 for (i = 0; i < nset; i++)
1305 {
1306 if (gmx_strcasecmp(leg[i], enm[set[i]].name) != 0)
1307 {
1308 fprintf(stderr, "\nWARNING energy set name mismatch %s!=%s\n",
1309 leg[i], enm[set[i]].name);
1310 }
1311 for (j = 0; j < nenergy; j++)
1312 {
1313 dE = eneset2[i][j] - edat->s[i].ener[j];
1314 sum += std::exp(-dE*beta);
1315 if (fp)
1316 {
1317 fprintf(fp, "%10g %10g %10g\n",
1318 time[j], dE, -BOLTZ*reftemp*std::log(sum/(j+1)) );
1319 }
1320 }
1321 aver = -BOLTZ*reftemp*std::log(sum/nenergy);
1322 fprintf(stdout, "%-24s %10g\n", leg[i], aver);
1323 }
1324 if (fp)
1325 {
1326 xvgrclose(fp);
1327 }
1328 sfree(fr);
1329 }
1330
1331
1332 static void do_dhdl(t_enxframe *fr, const t_inputrec *ir, FILE **fp_dhdl,
1333 const char *filename, gmx_bool bDp,
1334 int *blocks, int *hists, int *samples, int *nlambdas,
1335 const gmx_output_env_t *oenv)
1336 {
1337 const char *dhdl = "dH/d\\lambda", *deltag = "\\DeltaH", *lambda = "\\lambda";
1338 char title[STRLEN], label_x[STRLEN], label_y[STRLEN], legend[STRLEN];
1339 char buf[STRLEN];
1340 int nblock_hist = 0, nblock_dh = 0, nblock_dhcoll = 0;
1341 int i, j, k;
1342 /* coll data */
1343 double temp = 0, start_time = 0, delta_time = 0, start_lambda = 0;
1344 static int setnr = 0;
1345 double *native_lambda_vec = nullptr;
1346 const char **lambda_components = nullptr;
1347 int n_lambda_vec = 0;
1348 bool firstPass = true;
1349
1350 /* now count the blocks & handle the global dh data */
1351 for (i = 0; i < fr->nblock; i++)
1352 {
1353 if (fr->block[i].id == enxDHHIST)
1354 {
1355 nblock_hist++;
1356 }
1357 else if (fr->block[i].id == enxDH)
1358 {
1359 nblock_dh++;
1360 }
1361 else if (fr->block[i].id == enxDHCOLL)
1362 {
1363 nblock_dhcoll++;
1364 if ( (fr->block[i].nsub < 1) ||
1365 (fr->block[i].sub[0].type != xdr_datatype_double) ||
1366 (fr->block[i].sub[0].nr < 5))
1367 {
1368 gmx_fatal(FARGS, "Unexpected block data");
1369 }
1370
1371 /* read the data from the DHCOLL block */
1372 temp = fr->block[i].sub[0].dval[0];
1373 start_time = fr->block[i].sub[0].dval[1];
1374 delta_time = fr->block[i].sub[0].dval[2];
1375 start_lambda = fr->block[i].sub[0].dval[3];
1376 if (fr->block[i].nsub > 1)
1377 {
1378 if (firstPass)
1379 {
1380 n_lambda_vec = fr->block[i].sub[1].ival[1];
1381 snew(lambda_components, n_lambda_vec);
1382 snew(native_lambda_vec, n_lambda_vec);
1383 firstPass = false;
1384 }
1385 else
1386 {
1387 if (n_lambda_vec != fr->block[i].sub[1].ival[1])
1388 {
1389 gmx_fatal(FARGS,
1390 "Unexpected change of basis set in lambda");
1391 }
1392 }
1393 for (j = 0; j < n_lambda_vec; j++)
1394 {
1395 native_lambda_vec[j] = fr->block[i].sub[0].dval[5+j];
1396 lambda_components[j] =
1397 efpt_singular_names[fr->block[i].sub[1].ival[2+j]];
1398 }
1399 }
1400 }
1401 }
1402 // Clean up!
1403 sfree(native_lambda_vec);
1404 sfree(lambda_components);
1405 if (nblock_hist == 0 && nblock_dh == 0)
1406 {
1407 /* don't do anything */
1408 return;
1409 }
1410 if (nblock_hist > 0 && nblock_dh > 0)
1411 {
1412 gmx_fatal(FARGS, "This energy file contains both histogram dhdl data and non-histogram dhdl data. Don't know what to do.");
1413 }
1414 if (!*fp_dhdl)
1415 {
1416 if (nblock_dh > 0)
1417 {
1418 /* we have standard, non-histogram data --
1419 call open_dhdl to open the file */
1420 /* TODO this is an ugly hack that needs to be fixed: this will only
1421 work if the order of data is always the same and if we're
1422 only using the g_energy compiled with the mdrun that produced
1423 the ener.edr. */
1424 *fp_dhdl = open_dhdl(filename, ir, oenv);
1425 }
1426 else
1427 {
1428 sprintf(title, "N(%s)", deltag);
1429 sprintf(label_x, "%s (%s)", deltag, unit_energy);
1430 sprintf(label_y, "Samples");
1431 *fp_dhdl = xvgropen_type(filename, title, label_x, label_y, exvggtXNY, oenv);
1432 sprintf(buf, "T = %g (K), %s = %g", temp, lambda, start_lambda);
1433 xvgr_subtitle(*fp_dhdl, buf, oenv);
1434 }
1435 }
1436
1437 (*hists) += nblock_hist;
1438 (*blocks) += nblock_dh;
1439 (*nlambdas) = nblock_hist+nblock_dh;
1440
1441 /* write the data */
1442 if (nblock_hist > 0)
1443 {
1444 gmx_int64_t sum = 0;
1445 /* histograms */
1446 for (i = 0; i < fr->nblock; i++)
1447 {
1448 t_enxblock *blk = &(fr->block[i]);
1449 if (blk->id == enxDHHIST)
1450 {
1451 double foreign_lambda, dx;
1452 gmx_int64_t x0;
1453 int nhist, derivative;
1454
1455 /* check the block types etc. */
1456 if ( (blk->nsub < 2) ||
1457 (blk->sub[0].type != xdr_datatype_double) ||
1458 (blk->sub[1].type != xdr_datatype_int64) ||
1459 (blk->sub[0].nr < 2) ||
1460 (blk->sub[1].nr < 2) )
1461 {
1462 gmx_fatal(FARGS, "Unexpected block data in file");
1463 }
1464 foreign_lambda = blk->sub[0].dval[0];
1465 dx = blk->sub[0].dval[1];
1466 nhist = blk->sub[1].lval[0];
1467 derivative = blk->sub[1].lval[1];
1468 for (j = 0; j < nhist; j++)
1469 {
1470 const char *lg[1];
1471 x0 = blk->sub[1].lval[2+j];
1472
1473 if (!derivative)
1474 {
1475 sprintf(legend, "N(%s(%s=%g) | %s=%g)",
1476 deltag, lambda, foreign_lambda,
1477 lambda, start_lambda);
1478 }
1479 else
1480 {
1481 sprintf(legend, "N(%s | %s=%g)",
1482 dhdl, lambda, start_lambda);
1483 }
1484
1485 lg[0] = legend;
1486 xvgr_new_dataset(*fp_dhdl, setnr, 1, lg, oenv);
1487 setnr++;
1488 for (k = 0; k < blk->sub[j+2].nr; k++)
1489 {
1490 int hist;
1491 double xmin, xmax;
1492
1493 hist = blk->sub[j+2].ival[k];
1494 xmin = (x0+k)*dx;
1495 xmax = (x0+k+1)*dx;
1496 fprintf(*fp_dhdl, "%g %d\n%g %d\n", xmin, hist,
1497 xmax, hist);
1498 sum += hist;
1499 }
1500 /* multiple histogram data blocks in one histogram
1501 mean that the second one is the reverse of the first one:
1502 for dhdl derivatives, it's important to know both the
1503 maximum and minimum values */
1504 dx = -dx;
1505 }
1506 }
1507 }
1508 (*samples) += static_cast<int>(sum/nblock_hist);
1509 }
1510 else
1511 {
1512 /* raw dh */
1513 int len = 0;
1514
1515 for (i = 0; i < fr->nblock; i++)
1516 {
1517 t_enxblock *blk = &(fr->block[i]);
1518 if (blk->id == enxDH)
1519 {
1520 if (len == 0)
1521 {
1522 len = blk->sub[2].nr;
1523 }
1524 else
1525 {
1526 if (len != blk->sub[2].nr)
1527 {
1528 gmx_fatal(FARGS, "Length inconsistency in dhdl data");
1529 }
1530 }
1531 }
1532 }
1533 (*samples) += len;
1534
1535 for (i = 0; i < len; i++)
1536 {
1537 double time = start_time + delta_time*i;
1538
1539 fprintf(*fp_dhdl, "%.4f ", time);
1540
1541 for (j = 0; j < fr->nblock; j++)
1542 {
1543 t_enxblock *blk = &(fr->block[j]);
1544 if (blk->id == enxDH)
1545 {
1546 double value;
1547 if (blk->sub[2].type == xdr_datatype_float)
1548 {
1549 value = blk->sub[2].fval[i];
1550 }
1551 else
1552 {
1553 value = blk->sub[2].dval[i];
1554 }
1555 /* we need to decide which data type it is based on the count*/
1556
1557 if (j == 1 && ir->bExpanded)
1558 {
1559 fprintf(*fp_dhdl, "%4d", static_cast<int>(value)); /* if expanded ensembles and zero, this is a state value, it's an integer. We need a cleaner conditional than if j==1! */
1560 }
1561 else
1562 {
1563 if (bDp)
1564 {
1565 fprintf(*fp_dhdl, " %#.12g", value); /* print normal precision */
1566 }
1567 else
1568 {
1569 fprintf(*fp_dhdl, " %#.8g", value); /* print normal precision */
1570 }
1571 }
1572 }
1573 }
1574 fprintf(*fp_dhdl, "\n");
1575 }
1576 }
1577 }
1578
1579
1580 int gmx_energy(int argc, char *argv[])
1581 {
1582 const char *desc[] = {
1583 "[THISMODULE] extracts energy components",
1584 "from an energy file. The user is prompted to interactively",
1585 "select the desired energy terms.[PAR]",
1586
1587 "Average, RMSD, and drift are calculated with full precision from the",
1588 "simulation (see printed manual). Drift is calculated by performing",
1589 "a least-squares fit of the data to a straight line. The reported total drift",
1590 "is the difference of the fit at the first and last point.",
1591 "An error estimate of the average is given based on a block averages",
1592 "over 5 blocks using the full-precision averages. The error estimate",
1593 "can be performed over multiple block lengths with the options",
1594 "[TT]-nbmin[tt] and [TT]-nbmax[tt].",
1595 "[BB]Note[bb] that in most cases the energy files contains averages over all",
1596 "MD steps, or over many more points than the number of frames in",
1597 "energy file. This makes the [THISMODULE] statistics output more accurate",
1598 "than the [REF].xvg[ref] output. When exact averages are not present in the energy",
1599 "file, the statistics mentioned above are simply over the single, per-frame",
1600 "energy values.[PAR]",
1601
1602 "The term fluctuation gives the RMSD around the least-squares fit.[PAR]",
1603
1604 "Some fluctuation-dependent properties can be calculated provided",
1605 "the correct energy terms are selected, and that the command line option",
1606 "[TT]-fluct_props[tt] is given. The following properties",
1607 "will be computed:",
1608 "",
1609 "=============================== ===================",
1610 "Property Energy terms needed",
1611 "=============================== ===================",
1612 "Heat capacity C[SUB]p[sub] (NPT sims): Enthalpy, Temp",
1613 "Heat capacity C[SUB]v[sub] (NVT sims): Etot, Temp",
1614 "Thermal expansion coeff. (NPT): Enthalpy, Vol, Temp",
1615 "Isothermal compressibility: Vol, Temp",
1616 "Adiabatic bulk modulus: Vol, Temp",
1617 "=============================== ===================",
1618 "",
1619 "You always need to set the number of molecules [TT]-nmol[tt].",
1620 "The C[SUB]p[sub]/C[SUB]v[sub] computations do [BB]not[bb] include any corrections",
1621 "for quantum effects. Use the [gmx-dos] program if you need that (and you do).[PAR]"
1622
1623 "Option [TT]-odh[tt] extracts and plots the free energy data",
1624 "(Hamiltoian differences and/or the Hamiltonian derivative dhdl)",
1625 "from the [TT]ener.edr[tt] file.[PAR]",
1626
1627 "With [TT]-fee[tt] an estimate is calculated for the free-energy",
1628 "difference with an ideal gas state::",
1629 "",
1630 " [GRK]Delta[grk] A = A(N,V,T) - A[SUB]idealgas[sub](N,V,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln]",
1631 " [GRK]Delta[grk] G = G(N,p,T) - G[SUB]idealgas[sub](N,p,T) = kT [LN][CHEVRON][EXP]U[SUB]pot[sub]/kT[exp][chevron][ln]",
1632 "",
1633 "where k is Boltzmann's constant, T is set by [TT]-fetemp[tt] and",
1634 "the average is over the ensemble (or time in a trajectory).",
1635 "Note that this is in principle",
1636 "only correct when averaging over the whole (Boltzmann) ensemble",
1637 "and using the potential energy. This also allows for an entropy",
1638 "estimate using::",
1639 "",
1640 " [GRK]Delta[grk] S(N,V,T) = S(N,V,T) - S[SUB]idealgas[sub](N,V,T) = ([CHEVRON]U[SUB]pot[sub][chevron] - [GRK]Delta[grk] A)/T",
1641 " [GRK]Delta[grk] S(N,p,T) = S(N,p,T) - S[SUB]idealgas[sub](N,p,T) = ([CHEVRON]U[SUB]pot[sub][chevron] + pV - [GRK]Delta[grk] G)/T",
1642 "",
1643
1644 "When a second energy file is specified ([TT]-f2[tt]), a free energy",
1645 "difference is calculated::",
1646 "",
1647 " dF = -kT [LN][CHEVRON][EXP]-(E[SUB]B[sub]-E[SUB]A[sub])/kT[exp][chevron][SUB]A[sub][ln] ,",
1648 "",
1649 "where E[SUB]A[sub] and E[SUB]B[sub] are the energies from the first and second energy",
1650 "files, and the average is over the ensemble A. The running average",
1651 "of the free energy difference is printed to a file specified by [TT]-ravg[tt].",
1652 "[BB]Note[bb] that the energies must both be calculated from the same trajectory."
1653
1654 };
1655 static gmx_bool bSum = FALSE, bFee = FALSE, bPrAll = FALSE, bFluct = FALSE, bDriftCorr = FALSE;
1656 static gmx_bool bDp = FALSE, bMutot = FALSE, bOrinst = FALSE, bOvec = FALSE, bFluctProps = FALSE;
1657 static int nmol = 1, nbmin = 5, nbmax = 5;
1658 static real reftemp = 300.0, ezero = 0;
1659 t_pargs pa[] = {
1660 { "-fee", FALSE, etBOOL, {&bFee},
1661 "Do a free energy estimate" },
1662 { "-fetemp", FALSE, etREAL, {&reftemp},
1663 "Reference temperature for free energy calculation" },
1664 { "-zero", FALSE, etREAL, {&ezero},
1665 "Subtract a zero-point energy" },
1666 { "-sum", FALSE, etBOOL, {&bSum},
1667 "Sum the energy terms selected rather than display them all" },
1668 { "-dp", FALSE, etBOOL, {&bDp},
1669 "Print energies in high precision" },
1670 { "-nbmin", FALSE, etINT, {&nbmin},
1671 "Minimum number of blocks for error estimate" },
1672 { "-nbmax", FALSE, etINT, {&nbmax},
1673 "Maximum number of blocks for error estimate" },
1674 { "-mutot", FALSE, etBOOL, {&bMutot},
1675 "Compute the total dipole moment from the components" },
1676 { "-aver", FALSE, etBOOL, {&bPrAll},
1677 "Also print the exact average and rmsd stored in the energy frames (only when 1 term is requested)" },
1678 { "-nmol", FALSE, etINT, {&nmol},
1679 "Number of molecules in your sample: the energies are divided by this number" },
1680 { "-fluct_props", FALSE, etBOOL, {&bFluctProps},
1681 "Compute properties based on energy fluctuations, like heat capacity" },
1682 { "-driftcorr", FALSE, etBOOL, {&bDriftCorr},
1683 "Useful only for calculations of fluctuation properties. The drift in the observables will be subtracted before computing the fluctuation properties."},
1684 { "-fluc", FALSE, etBOOL, {&bFluct},
1685 "Calculate autocorrelation of energy fluctuations rather than energy itself" },
1686 { "-orinst", FALSE, etBOOL, {&bOrinst},
1687 "Analyse instantaneous orientation data" },
1688 { "-ovec", FALSE, etBOOL, {&bOvec},
1689 "Also plot the eigenvectors with [TT]-oten[tt]" }
1690 };
1691 static const char *setnm[] = {
1692 "Pres-XX", "Pres-XY", "Pres-XZ", "Pres-YX", "Pres-YY",
1693 "Pres-YZ", "Pres-ZX", "Pres-ZY", "Pres-ZZ", "Temperature",
1694 "Volume", "Pressure"
1695 };
1696
1697 FILE *out = nullptr;
1698 FILE *fp_dhdl = nullptr;
1699 ener_file_t fp;
1700 int timecheck = 0;
1701 enerdata_t edat;
1702 gmx_enxnm_t *enm = nullptr;
1703 t_enxframe *frame, *fr = nullptr;
1704 int cur = 0;
1705 #define NEXT (1-cur)
1706 int nre, nfr;
1707 gmx_int64_t start_step;
1708 real start_t;
1709 gmx_bool bDHDL;
1710 gmx_bool bFoundStart, bCont, bVisco;
1711 double sum, dbl;
1712 double *time = nullptr;
1713 real Vaver;
1714 int *set = nullptr, i, j, nset, sss;
1715 gmx_bool *bIsEner = nullptr;
1716 char **leg = nullptr;
1717 char buf[256];
1718 gmx_output_env_t *oenv;
1719 int dh_blocks = 0, dh_hists = 0, dh_samples = 0, dh_lambdas = 0;
1720
1721 t_filenm fnm[] = {
1722 { efEDR, "-f", nullptr, ffREAD },
1723 { efEDR, "-f2", nullptr, ffOPTRD },
1724 { efTPR, "-s", nullptr, ffOPTRD },
1725 { efXVG, "-o", "energy", ffWRITE },
1726 { efXVG, "-viol", "violaver", ffOPTWR },
1727 { efXVG, "-pairs", "pairs", ffOPTWR },
1728 { efXVG, "-corr", "enecorr", ffOPTWR },
1729 { efXVG, "-vis", "visco", ffOPTWR },
1730 { efXVG, "-evisco", "evisco", ffOPTWR },
1731 { efXVG, "-eviscoi", "eviscoi", ffOPTWR },
1732 { efXVG, "-ravg", "runavgdf", ffOPTWR },
1733 { efXVG, "-odh", "dhdl", ffOPTWR }
1734 };
1735 #define NFILE asize(fnm)
1736 int npargs;
1737 t_pargs *ppa;
1738
1739 npargs = asize(pa);
1740 ppa = add_acf_pargs(&npargs, pa);
1741 if (!parse_common_args(&argc, argv,
1742 PCA_CAN_VIEW | PCA_CAN_BEGIN | PCA_CAN_END,
1743 NFILE, fnm, npargs, ppa, asize(desc), desc, 0, nullptr, &oenv))
1744 {
1745 sfree(ppa);
1746 return 0;
1747 }
1748
1749 bDHDL = opt2bSet("-odh", NFILE, fnm);
1750
1751 nset = 0;
1752
1753 snew(frame, 2);
1754 fp = open_enx(ftp2fn(efEDR, NFILE, fnm), "r");
1755 do_enxnms(fp, &nre, &enm);
1756
1757 Vaver = -1;
1758
1759 bVisco = opt2bSet("-vis", NFILE, fnm);
1760
1761 t_inputrec irInstance;
1762 t_inputrec *ir = &irInstance;
1763
1764 if (!bDHDL)
1765 {
1766 if (bVisco)
1767 {
1768 nset = asize(setnm);
1769 snew(set, nset);
1770 /* This is nasty code... To extract Pres tensor, Volume and Temperature */
1771 for (j = 0; j < nset; j++)
1772 {
1773 for (i = 0; i < nre; i++)
1774 {
1775 if (std::strstr(enm[i].name, setnm[j]))
1776 {
1777 set[j] = i;
1778 break;
1779 }
1780 }
1781 if (i == nre)
1782 {
1783 if (gmx_strcasecmp(setnm[j], "Volume") == 0)
1784 {
1785 printf("Enter the box volume (" unit_volume "): ");
1786 if (1 != scanf("%lf", &dbl))
1787 {
1788 gmx_fatal(FARGS, "Error reading user input");
1789 }
1790 Vaver = dbl;
1791 }
1792 else
1793 {
1794 gmx_fatal(FARGS, "Could not find term %s for viscosity calculation",
1795 setnm[j]);
1796 }
1797 }
1798 }
1799 }
1800 else
1801 {
1802 set = select_by_name(nre, enm, &nset);
1803 }
1804 /* Print all the different units once */
1805 sprintf(buf, "(%s)", enm[set[0]].unit);
1806 for (i = 1; i < nset; i++)
1807 {
1808 for (j = 0; j < i; j++)
1809 {
1810 if (std::strcmp(enm[set[i]].unit, enm[set[j]].unit) == 0)
1811 {
1812 break;
1813 }
1814 }
1815 if (j == i)
1816 {
1817 std::strcat(buf, ", (");
1818 std::strcat(buf, enm[set[i]].unit);
1819 std::strcat(buf, ")");
1820 }
1821 }
1822 out = xvgropen(opt2fn("-o", NFILE, fnm), "GROMACS Energies", "Time (ps)", buf,
1823 oenv);
1824
1825 snew(leg, nset+1);
1826 for (i = 0; (i < nset); i++)
1827 {
1828 leg[i] = enm[set[i]].name;
1829 }
1830 if (bSum)
1831 {
1832 leg[nset] = gmx_strdup("Sum");
1833 xvgr_legend(out, nset+1, (const char**)leg, oenv);
1834 }
1835 else
1836 {
1837 xvgr_legend(out, nset, (const char**)leg, oenv);
1838 }
1839
1840 snew(bIsEner, nset);
1841 for (i = 0; (i < nset); i++)
1842 {
1843 bIsEner[i] = FALSE;
1844 for (j = 0; (j <= F_ETOT); j++)
1845 {
1846 bIsEner[i] = bIsEner[i] ||
1847 (gmx_strcasecmp(interaction_function[j].longname, leg[i]) == 0);
1848 }
1849 }
1850 if (bPrAll && nset > 1)
1851 {
1852 gmx_fatal(FARGS, "Printing averages can only be done when a single set is selected");
1853 }
1854
1855 }
1856 else if (bDHDL)
1857 {
1858 get_dhdl_parms(ftp2fn(efTPR, NFILE, fnm), ir);
1859 }
1860
1861 /* Initiate energies and set them to zero */
1862 edat.nsteps = 0;
1863 edat.npoints = 0;
1864 edat.nframes = 0;
1865 edat.step = nullptr;
1866 edat.steps = nullptr;
1867 edat.points = nullptr;
1868 edat.bHaveSums = TRUE;
1869 snew(edat.s, nset);
1870
1871 /* Initiate counters */
1872 bFoundStart = FALSE;
1873 start_step = 0;
1874 start_t = 0;
1875 do
1876 {
1877 /* This loop searches for the first frame (when -b option is given),
1878 * or when this has been found it reads just one energy frame
1879 */
1880 do
1881 {
1882 bCont = do_enx(fp, &(frame[NEXT]));
1883 if (bCont)
1884 {
1885 timecheck = check_times(frame[NEXT].t);
1886 }
1887 }
1888 while (bCont && (timecheck < 0));
1889
1890 if ((timecheck == 0) && bCont)
1891 {
1892 /* We read a valid frame, so we can use it */
1893 fr = &(frame[NEXT]);
1894
1895 if (fr->nre > 0)
1896 {
1897 /* The frame contains energies, so update cur */
1898 cur = NEXT;
1899
1900 if (edat.nframes % 1000 == 0)
1901 {
1902 srenew(edat.step, edat.nframes+1000);
1903 std::memset(&(edat.step[edat.nframes]), 0, 1000*sizeof(edat.step[0]));
1904 srenew(edat.steps, edat.nframes+1000);
1905 std::memset(&(edat.steps[edat.nframes]), 0, 1000*sizeof(edat.steps[0]));
1906 srenew(edat.points, edat.nframes+1000);
1907 std::memset(&(edat.points[edat.nframes]), 0, 1000*sizeof(edat.points[0]));
1908
1909 for (i = 0; i < nset; i++)
1910 {
1911 srenew(edat.s[i].ener, edat.nframes+1000);
1912 std::memset(&(edat.s[i].ener[edat.nframes]), 0,
1913 1000*sizeof(edat.s[i].ener[0]));
1914 srenew(edat.s[i].es, edat.nframes+1000);
1915 std::memset(&(edat.s[i].es[edat.nframes]), 0,
1916 1000*sizeof(edat.s[i].es[0]));
1917 }
1918 }
1919
1920 nfr = edat.nframes;
1921 edat.step[nfr] = fr->step;
1922
1923 if (!bFoundStart)
1924 {
1925 bFoundStart = TRUE;
1926 /* Initiate the previous step data */
1927 start_step = fr->step;
1928 start_t = fr->t;
1929 /* Initiate the energy sums */
1930 edat.steps[nfr] = 1;
1931 edat.points[nfr] = 1;
1932 for (i = 0; i < nset; i++)
1933 {
1934 sss = set[i];
1935 edat.s[i].es[nfr].sum = fr->ener[sss].e;
1936 edat.s[i].es[nfr].sum2 = 0;
1937 }
1938 edat.nsteps = 1;
1939 edat.npoints = 1;
1940 }
1941 else
1942 {
1943 edat.steps[nfr] = fr->nsteps;
1944
1945 if (fr->nsum <= 1)
1946 {
1947 /* mdrun only calculated the energy at energy output
1948 * steps. We don't need to check step intervals.
1949 */
1950 edat.points[nfr] = 1;
1951 for (i = 0; i < nset; i++)
1952 {
1953 sss = set[i];
1954 edat.s[i].es[nfr].sum = fr->ener[sss].e;
1955 edat.s[i].es[nfr].sum2 = 0;
1956 }
1957 edat.npoints += 1;
1958 edat.bHaveSums = FALSE;
1959 }
1960 else if (fr->step - start_step + 1 == edat.nsteps + fr->nsteps)
1961 {
1962 /* We have statistics to the previous frame */
1963 edat.points[nfr] = fr->nsum;
1964 for (i = 0; i < nset; i++)
1965 {
1966 sss = set[i];
1967 edat.s[i].es[nfr].sum = fr->ener[sss].esum;
1968 edat.s[i].es[nfr].sum2 = fr->ener[sss].eav;
1969 }
1970 edat.npoints += fr->nsum;
1971 }
1972 else
1973 {
1974 /* The interval does not match fr->nsteps:
1975 * can not do exact averages.
1976 */
1977 edat.bHaveSums = FALSE;
1978 }
1979
1980 edat.nsteps = fr->step - start_step + 1;
1981 }
1982 for (i = 0; i < nset; i++)
1983 {
1984 edat.s[i].ener[nfr] = fr->ener[set[i]].e;
1985 }
1986 }
1987 /*
1988 * Store energies for analysis afterwards...
1989 */
1990 if (!bDHDL && (fr->nre > 0))
1991 {
1992 if (edat.nframes % 1000 == 0)
1993 {
1994 srenew(time, edat.nframes+1000);
1995 }
1996 time[edat.nframes] = fr->t;
1997 edat.nframes++;
1998 }
1999 if (bDHDL)
2000 {
2001 do_dhdl(fr, ir, &fp_dhdl, opt2fn("-odh", NFILE, fnm), bDp, &dh_blocks, &dh_hists, &dh_samples, &dh_lambdas, oenv);
2002 }
2003
2004 /*******************************************
2005 * E N E R G I E S
2006 *******************************************/
2007 else
2008 {
2009 if (fr->nre > 0)
2010 {
2011 if (bPrAll)
2012 {
2013 /* We skip frames with single points (usually only the first frame),
2014 * since they would result in an average plot with outliers.
2015 */
2016 if (fr->nsum > 1)
2017 {
2018 print_time(out, fr->t);
2019 print1(out, bDp, fr->ener[set[0]].e);
2020 print1(out, bDp, fr->ener[set[0]].esum/fr->nsum);
2021 print1(out, bDp, std::sqrt(fr->ener[set[0]].eav/fr->nsum));
2022 fprintf(out, "\n");
2023 }
2024 }
2025 else
2026 {
2027 print_time(out, fr->t);
2028 if (bSum)
2029 {
2030 sum = 0;
2031 for (i = 0; i < nset; i++)
2032 {
2033 sum += fr->ener[set[i]].e;
2034 }
2035 print1(out, bDp, sum/nmol-ezero);
2036 }
2037 else
2038 {
2039 for (i = 0; (i < nset); i++)
2040 {
2041 if (bIsEner[i])
2042 {
2043 print1(out, bDp, (fr->ener[set[i]].e)/nmol-ezero);
2044 }
2045 else
2046 {
2047 print1(out, bDp, fr->ener[set[i]].e);
2048 }
2049 }
2050 }
2051 fprintf(out, "\n");
2052 }
2053 }
2054 }
2055 }
2056 }
2057 while (bCont && (timecheck == 0));
2058
2059 fprintf(stderr, "\n");
2060 done_ener_file(fp);
2061 if (out)
2062 {
2063 xvgrclose(out);
2064 }
2065
2066 if (bDHDL)
2067 {
2068 if (fp_dhdl)
2069 {
2070 gmx_fio_fclose(fp_dhdl);
2071 printf("\n\nWrote %d lambda values with %d samples as ",
2072 dh_lambdas, dh_samples);
2073 if (dh_hists > 0)
2074 {
2075 printf("%d dH histograms ", dh_hists);
2076 }
2077 if (dh_blocks > 0)
2078 {
2079 printf("%d dH data blocks ", dh_blocks);
2080 }
2081 printf("to %s\n", opt2fn("-odh", NFILE, fnm));
2082
2083 }
2084 else
2085 {
2086 gmx_fatal(FARGS, "No dH data in %s\n", opt2fn("-f", NFILE, fnm));
2087 }
2088
2089 }
2090 else
2091 {
2092 double dt = (frame[cur].t-start_t)/(edat.nframes-1);
2093 analyse_ener(opt2bSet("-corr", NFILE, fnm), opt2fn("-corr", NFILE, fnm),
2094 opt2fn("-evisco", NFILE, fnm), opt2fn("-eviscoi", NFILE, fnm),
2095 bFee, bSum, bFluct,
2096 bVisco, opt2fn("-vis", NFILE, fnm),
2097 nmol,
2098 start_step, start_t, frame[cur].step, frame[cur].t,
2099 reftemp, &edat,
2100 nset, set, bIsEner, leg, enm, Vaver, ezero, nbmin, nbmax,
2101 oenv);
2102 if (bFluctProps)
2103 {
2104 calc_fluctuation_props(stdout, bDriftCorr, dt, nset, nmol, leg, &edat,
2105 nbmin, nbmax);
2106 }
2107 }
2108 if (opt2bSet("-f2", NFILE, fnm))
2109 {
2110 fec(opt2fn("-f2", NFILE, fnm), opt2fn("-ravg", NFILE, fnm),
2111 reftemp, nset, set, leg, &edat, time, oenv);
2112 }
2113 // Clean up!
2114 done_enerdata_t(nset, &edat);
2115 sfree(time);
2116 free_enxframe(&frame[0]);
2117 free_enxframe(&frame[1]);
2118 sfree(frame);
2119 free_enxnms(nre, enm);
2120 sfree(ppa);
2121 sfree(set);
2122 sfree(leg);
2123 sfree(bIsEner);
2124 {
2125 const char *nxy = "-nxy";
2126
2127 do_view(oenv, opt2fn("-o", NFILE, fnm), nxy);
2128 do_view(oenv, opt2fn_null("-ravg", NFILE, fnm), nxy);
2129 do_view(oenv, opt2fn_null("-odh", NFILE, fnm), nxy);
2130 }
2131 output_env_done(oenv);
2132
2133 return 0;
2134 }
The ultimate molecular dynamics simulation package