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