/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx

Bug Summary

File:	gromacs/gmxana/gmx_energy.c
Location:	line 1038, column 5
Description:	Value stored to 'et' is never read

Annotated Source Code

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