Bug Summary

File:gromacs/gmxana/gmx_eneconv.c
Location:line 586, column 5
Description:Value stored to 'nfile' 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/smalloc.h"
47#include "gromacs/commandline/pargs.h"
48#include "disre.h"
49#include "names.h"
50#include "macros.h"
51#include "gromacs/utility/fatalerror.h"
52#include "gromacs/fileio/enxio.h"
53#include "gromacs/math/vec.h"
54#include "gmx_ana.h"
55#include "gromacs/fileio/trxio.h"
56
57#define TIME_EXPLICIT0 0
58#define TIME_CONTINUE1 1
59#define TIME_LAST2 2
60#ifndef FLT_MAX1e36
61#define FLT_MAX1e36 1e36
62#endif
63
64static int *select_it(int nre, gmx_enxnm_t *nm, int *nset)
65{
66 gmx_bool *bE;
67 int n, k, j, i;
68 int *set;
69 gmx_bool bVerbose = TRUE1;
70
71 if ((getenv("VERBOSE")) != NULL((void*)0))
72 {
73 bVerbose = FALSE0;
74 }
75
76 fprintf(stderrstderr, "Select the terms you want to scale from the following list\n");
77 fprintf(stderrstderr, "End your selection with 0\n");
78
79 if (bVerbose)
80 {
81 for (k = 0; (k < nre); )
82 {
83 for (j = 0; (j < 4) && (k < nre); j++, k++)
84 {
85 fprintf(stderrstderr, " %3d=%14s", k+1, nm[k].name);
86 }
87 fprintf(stderrstderr, "\n");
88 }
89 }
90
91 snew(bE, nre)(bE) = save_calloc("bE", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 91, (nre), sizeof(*(bE)));
92 do
93 {
94 if (1 != scanf("%d", &n))
95 {
96 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 96, "Cannot read energy term");
97 }
98 if ((n > 0) && (n <= nre))
99 {
100 bE[n-1] = TRUE1;
101 }
102 }
103 while (n != 0);
104
105 snew(set, nre)(set) = save_calloc("set", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 105, (nre), sizeof(*(set)));
106 for (i = (*nset) = 0; (i < nre); i++)
107 {
108 if (bE[i])
109 {
110 set[(*nset)++] = i;
111 }
112 }
113
114 sfree(bE)save_free("bE", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 114, (bE));
115
116 return set;
117}
118
119static void sort_files(char **fnms, real *settime, int nfile)
120{
121 int i, j, minidx;
122 real timeswap;
123 char *chptr;
124
125 for (i = 0; i < nfile; i++)
126 {
127 minidx = i;
128 for (j = i+1; j < nfile; j++)
129 {
130 if (settime[j] < settime[minidx])
131 {
132 minidx = j;
133 }
134 }
135 if (minidx != i)
136 {
137 timeswap = settime[i];
138 settime[i] = settime[minidx];
139 settime[minidx] = timeswap;
140 chptr = fnms[i];
141 fnms[i] = fnms[minidx];
142 fnms[minidx] = chptr;
143 }
144 }
145}
146
147
148static int scan_ene_files(char **fnms, int nfiles,
149 real *readtime, real *timestep, int *nremax)
150{
151 /* Check number of energy terms and start time of all files */
152 int f, i, nre, nremin = 0, nresav = 0;
153 ener_file_t in;
154 real t1, t2;
155 char inputstring[STRLEN4096];
156 gmx_enxnm_t *enm;
157 t_enxframe *fr;
158
159 snew(fr, 1)(fr) = save_calloc("fr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 159, (1), sizeof(*(fr)));
160
161 for (f = 0; f < nfiles; f++)
162 {
163 in = open_enx(fnms[f], "r");
164 enm = NULL((void*)0);
165 do_enxnms(in, &nre, &enm);
166
167 if (f == 0)
168 {
169 nresav = nre;
170 nremin = nre;
171 *nremax = nre;
172 do_enx(in, fr);
173 t1 = fr->t;
174 do_enx(in, fr);
175 t2 = fr->t;
176 *timestep = t2-t1;
177 readtime[f] = t1;
178 close_enx(in);
179 }
180 else
181 {
182 nremin = min(nremin, fr->nre)(((nremin) < (fr->nre)) ? (nremin) : (fr->nre) );
183 *nremax = max(*nremax, fr->nre)(((*nremax) > (fr->nre)) ? (*nremax) : (fr->nre) );
184 if (nre != nresav)
185 {
186 fprintf(stderrstderr,
187 "Energy files don't match, different number of energies:\n"
188 " %s: %d\n %s: %d\n", fnms[f-1], nresav, fnms[f], fr->nre);
189 fprintf(stderrstderr,
190 "\nContinue conversion using only the first %d terms (n/y)?\n"
191 "(you should be sure that the energy terms match)\n", nremin);
192 if (NULL((void*)0) == fgets(inputstring, STRLEN4096-1, stdinstdin))
193 {
194 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 194, "Error reading user input");
195 }
196 if (inputstring[0] != 'y' && inputstring[0] != 'Y')
197 {
198 fprintf(stderrstderr, "Will not convert\n");
199 exit(0);
200 }
201 nresav = fr->nre;
202 }
203 do_enx(in, fr);
204 readtime[f] = fr->t;
205 close_enx(in);
206 }
207 fprintf(stderrstderr, "\n");
208 free_enxnms(nre, enm);
209 }
210
211 free_enxframe(fr);
212 sfree(fr)save_free("fr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 212, (fr));
213
214 return nremin;
215}
216
217
218static void edit_files(char **fnms, int nfiles, real *readtime,
219 real *settime, int *cont_type, gmx_bool bSetTime, gmx_bool bSort)
220{
221 int i;
222 gmx_bool ok;
223 char inputstring[STRLEN4096], *chptr;
224
225 if (bSetTime)
226 {
227 if (nfiles == 1)
228 {
229 fprintf(stderrstderr, "\n\nEnter the new start time:\n\n");
230 }
231 else
232 {
233 fprintf(stderrstderr, "\n\nEnter the new start time for each file.\n"
234 "There are two special options, both disables sorting:\n\n"
235 "c (continue) - The start time is taken from the end\n"
236 "of the previous file. Use it when your continuation run\n"
237 "restarts with t=0 and there is no overlap.\n\n"
238 "l (last) - The time in this file will be changed the\n"
239 "same amount as in the previous. Use it when the time in the\n"
240 "new run continues from the end of the previous one,\n"
241 "since this takes possible overlap into account.\n\n");
242 }
243
244 fprintf(stderrstderr, " File Current start New start\n"
245 "---------------------------------------------------------\n");
246
247 for (i = 0; i < nfiles; i++)
248 {
249 fprintf(stderrstderr, "%25s %10.3f ", fnms[i], readtime[i]);
250 ok = FALSE0;
251 do
252 {
253 if (NULL((void*)0) == fgets(inputstring, STRLEN4096-1, stdinstdin))
254 {
255 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 255, "Error reading user input");
256 }
257 inputstring[strlen(inputstring)-1] = 0;
258
259 if (inputstring[0] == 'c' || inputstring[0] == 'C')
260 {
261 cont_type[i] = TIME_CONTINUE1;
262 bSort = FALSE0;
263 ok = TRUE1;
264 settime[i] = FLT_MAX1e36;
265 }
266 else if (inputstring[0] == 'l' ||
267 inputstring[0] == 'L')
268 {
269 cont_type[i] = TIME_LAST2;
270 bSort = FALSE0;
271 ok = TRUE1;
272 settime[i] = FLT_MAX1e36;
273 }
274 else
275 {
276 settime[i] = strtod(inputstring, &chptr);
277 if (chptr == inputstring)
278 {
279 fprintf(stderrstderr, "Try that again: ");
280 }
281 else
282 {
283 cont_type[i] = TIME_EXPLICIT0;
284 ok = TRUE1;
285 }
286 }
287 }
288 while (!ok);
289 }
290 if (cont_type[0] != TIME_EXPLICIT0)
291 {
292 cont_type[0] = TIME_EXPLICIT0;
293 settime[0] = 0;
294 }
295 }
296 else
297 {
298 for (i = 0; i < nfiles; i++)
299 {
300 settime[i] = readtime[i];
301 }
302 }
303
304 if (bSort && (nfiles > 1))
305 {
306 sort_files(fnms, settime, nfiles);
307 }
308 else
309 {
310 fprintf(stderrstderr, "Sorting disabled.\n");
311 }
312
313
314 /* Write out the new order and start times */
315 fprintf(stderrstderr, "\nSummary of files and start times used:\n\n"
316 " File Start time\n"
317 "-----------------------------------------\n");
318 for (i = 0; i < nfiles; i++)
319 {
320 switch (cont_type[i])
321 {
322 case TIME_EXPLICIT0:
323 fprintf(stderrstderr, "%25s %10.3f\n", fnms[i], settime[i]);
324 break;
325 case TIME_CONTINUE1:
326 fprintf(stderrstderr, "%25s Continue from end of last file\n", fnms[i]);
327 break;
328 case TIME_LAST2:
329 fprintf(stderrstderr, "%25s Change by same amount as last file\n", fnms[i]);
330 break;
331 }
332 }
333 fprintf(stderrstderr, "\n");
334
335 settime[nfiles] = FLT_MAX1e36;
336 cont_type[nfiles] = TIME_EXPLICIT0;
337 readtime[nfiles] = FLT_MAX1e36;
338}
339
340
341static void copy_ee(t_energy *src, t_energy *dst, int nre)
342{
343 int i;
344
345 for (i = 0; i < nre; i++)
346 {
347 dst[i].e = src[i].e;
348 dst[i].esum = src[i].esum;
349 dst[i].eav = src[i].eav;
350 }
351}
352
353static void update_ee(t_energy *lastee, gmx_int64_t laststep,
354 t_energy *startee, gmx_int64_t startstep,
355 t_energy *ee, int step,
356 t_energy *outee, int nre)
357{
358 int i;
359 double sigmacorr, nom, denom;
360 double prestart_esum;
361 double prestart_sigma;
362
363 for (i = 0; i < nre; i++)
364 {
365 outee[i].e = ee[i].e;
366 /* add statistics from earlier file if present */
367 if (laststep > 0)
368 {
369 outee[i].esum = lastee[i].esum+ee[i].esum;
370 nom = (lastee[i].esum*(step+1)-ee[i].esum*(laststep));
371 denom = ((step+1.0)*(laststep)*(step+1.0+laststep));
372 sigmacorr = nom*nom/denom;
373 outee[i].eav = lastee[i].eav+ee[i].eav+sigmacorr;
374 }
375 else
376 {
377 /* otherwise just copy to output */
378 outee[i].esum = ee[i].esum;
379 outee[i].eav = ee[i].eav;
380 }
381
382 /* if we didnt start to write at the first frame
383 * we must compensate the statistics for this
384 * there are laststep frames in the earlier file
385 * and step+1 frames in this one.
386 */
387 if (startstep > 0)
388 {
389 gmx_int64_t q = laststep+step;
390 gmx_int64_t p = startstep+1;
391 prestart_esum = startee[i].esum-startee[i].e;
392 sigmacorr = prestart_esum-(p-1)*startee[i].e;
393 prestart_sigma = startee[i].eav-
394 sigmacorr*sigmacorr/(p*(p-1));
395 sigmacorr = prestart_esum/(p-1)-
396 outee[i].esum/(q);
397 outee[i].esum -= prestart_esum;
398 if (q-p+1 > 0)
399 {
400 outee[i].eav = outee[i].eav-prestart_sigma-
401 sigmacorr*sigmacorr*((p-1)*q)/(q-p+1);
402 }
403 }
404
405 if ((outee[i].eav/(laststep+step+1)) < (GMX_REAL_EPS5.96046448E-08))
406 {
407 outee[i].eav = 0;
408 }
409 }
410}
411
412static void update_ee_sum(int nre,
413 gmx_int64_t *ee_sum_step,
414 gmx_int64_t *ee_sum_nsteps,
415 gmx_int64_t *ee_sum_nsum,
416 t_energy *ee_sum,
417 t_enxframe *fr, int out_step)
418{
419 gmx_int64_t nsteps, nsum, fr_nsum;
420 int i;
421
422 nsteps = *ee_sum_nsteps;
423 nsum = *ee_sum_nsum;
424
425 fr_nsum = fr->nsum;
426 if (fr_nsum == 0)
427 {
428 fr_nsum = 1;
429 }
430
431 if (nsteps == 0)
432 {
433 if (fr_nsum == 1)
434 {
435 for (i = 0; i < nre; i++)
436 {
437 ee_sum[i].esum = fr->ener[i].e;
438 ee_sum[i].eav = 0;
439 }
440 }
441 else
442 {
443 for (i = 0; i < nre; i++)
444 {
445 ee_sum[i].esum = fr->ener[i].esum;
446 ee_sum[i].eav = fr->ener[i].eav;
447 }
448 }
449 nsteps = fr->nsteps;
450 nsum = fr_nsum;
451 }
452 else if (out_step + *ee_sum_nsum - *ee_sum_step == nsteps + fr->nsteps)
453 {
454 if (fr_nsum == 1)
455 {
456 for (i = 0; i < nre; i++)
457 {
458 ee_sum[i].eav +=
459 dsqr(ee_sum[i].esum/nsum
460 - (ee_sum[i].esum + fr->ener[i].e)/(nsum + 1))*nsum*(nsum + 1);
461 ee_sum[i].esum += fr->ener[i].e;
462 }
463 }
464 else
465 {
466 for (i = 0; i < fr->nre; i++)
467 {
468 ee_sum[i].eav +=
469 fr->ener[i].eav +
470 dsqr(ee_sum[i].esum/nsum
471 - (ee_sum[i].esum + fr->ener[i].esum)/(nsum + fr->nsum))*
472 nsum*(nsum + fr->nsum)/(double)fr->nsum;
473 ee_sum[i].esum += fr->ener[i].esum;
474 }
475 }
476 nsteps += fr->nsteps;
477 nsum += fr_nsum;
478 }
479 else
480 {
481 if (fr->nsum != 0)
482 {
483 fprintf(stderrstderr, "\nWARNING: missing energy sums at time %f\n", fr->t);
484 }
485 nsteps = 0;
486 nsum = 0;
487 }
488
489 *ee_sum_step = out_step;
490 *ee_sum_nsteps = nsteps;
491 *ee_sum_nsum = nsum;
492}
493
494int gmx_eneconv(int argc, char *argv[])
495{
496 const char *desc[] = {
497 "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[BR]",
498 "Concatenates several energy files in sorted order.",
499 "In the case of double time frames, the one",
500 "in the later file is used. By specifying [TT]-settime[tt] you will be",
501 "asked for the start time of each file. The input files are taken",
502 "from the command line,",
503 "such that the command [TT]gmx eneconv -f *.edr -o fixed.edr[tt] should do",
504 "the trick. [PAR]",
505 "With [IT]one file[it] specified for [TT]-f[tt]:[BR]",
506 "Reads one energy file and writes another, applying the [TT]-dt[tt],",
507 "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
508 "converting to a different format if necessary (indicated by file",
509 "extentions).[PAR]",
510 "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
511 "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
512 };
513 const char *bugs[] = {
514 "When combining trajectories the sigma and E^2 (necessary for statistics) are not updated correctly. Only the actual energy is correct. One thus has to compute statistics in another way."
515 };
516 ener_file_t in = NULL((void*)0), out = NULL((void*)0);
517 gmx_enxnm_t *enm = NULL((void*)0);
518#if 0
519 ener_file_t in, out = NULL((void*)0);
520 gmx_enxnm_t *enm = NULL((void*)0);
521#endif
522 t_enxframe *fr, *fro;
523 gmx_int64_t ee_sum_step = 0, ee_sum_nsteps, ee_sum_nsum;
524 t_energy *ee_sum;
525 gmx_int64_t lastfilestep, laststep, startstep, startstep_file = 0;
526 int noutfr;
527 int nre, nremax, this_nre, nfile, f, i, j, kkk, nset, *set = NULL((void*)0);
528 double last_t;
529 char **fnms;
530 real *readtime, *settime, timestep, t1, tadjust;
531 char inputstring[STRLEN4096], *chptr, buf[22], buf2[22], buf3[22];
532 gmx_bool ok;
533 int *cont_type;
534 gmx_bool bNewFile, bFirst, bNewOutput;
535 output_env_t oenv;
536 gmx_bool warned_about_dh = FALSE0;
537 t_enxblock *blocks = NULL((void*)0);
538 int nblocks = 0;
539 int nblocks_alloc = 0;
540
541 t_filenm fnm[] = {
542 { efEDR, "-f", NULL((void*)0), ffRDMULT(1<<1 | 1<<5) },
543 { efEDR, "-o", "fixed", ffWRITE1<<2 },
544 };
545
546#define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0])))
547 gmx_bool bWrite;
548 static real delta_t = 0.0, toffset = 0, scalefac = 1;
549 static gmx_bool bSetTime = FALSE0;
550 static gmx_bool bSort = TRUE1, bError = TRUE1;
551 static real begin = -1;
552 static real end = -1;
553 gmx_bool remove_dh = FALSE0;
554
555 t_pargs pa[] = {
556 { "-b", FALSE0, etREAL, {&begin},
557 "First time to use"},
558 { "-e", FALSE0, etREAL, {&end},
559 "Last time to use"},
560 { "-dt", FALSE0, etREAL, {&delta_t},
561 "Only write out frame when t MOD dt = offset" },
562 { "-offset", FALSE0, etREAL, {&toffset},
563 "Time offset for [TT]-dt[tt] option" },
564 { "-settime", FALSE0, etBOOL, {&bSetTime},
565 "Change starting time interactively" },
566 { "-sort", FALSE0, etBOOL, {&bSort},
567 "Sort energy files (not frames)"},
568 { "-rmdh", FALSE0, etBOOL, {&remove_dh},
569 "Remove free energy block data" },
570 { "-scalefac", FALSE0, etREAL, {&scalefac},
571 "Multiply energy component by this factor" },
572 { "-error", FALSE0, etBOOL, {&bError},
573 "Stop on errors in the file" }
574 };
575
576 if (!parse_common_args(&argc, argv, PCA_BE_NICE(1<<13), NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))),
577 pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), desc, asize(bugs)((int)(sizeof(bugs)/sizeof((bugs)[0]))), bugs, &oenv))
578 {
579 return 0;
580 }
581 tadjust = 0;
582 nremax = 0;
583 nset = 0;
584 timestep = 0.0;
585 snew(fnms, argc)(fnms) = save_calloc("fnms", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 585, (argc), sizeof(*(fnms)));
586 nfile = 0;
Value stored to 'nfile' is never read
587 lastfilestep = 0;
588 laststep = startstep = 0;
589
590 nfile = opt2fns(&fnms, "-f", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm);
591
592 if (!nfile)
593 {
594 gmx_fatal(FARGS0, "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 594, "No input files!");
595 }
596
597 snew(settime, nfile+1)(settime) = save_calloc("settime", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 597, (nfile+1), sizeof(*(settime)));
598 snew(readtime, nfile+1)(readtime) = save_calloc("readtime", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 598, (nfile+1), sizeof(*(readtime)));
599 snew(cont_type, nfile+1)(cont_type) = save_calloc("cont_type", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 599, (nfile+1), sizeof(*(cont_type)));
600
601 nre = scan_ene_files(fnms, nfile, readtime, &timestep, &nremax);
602 edit_files(fnms, nfile, readtime, settime, cont_type, bSetTime, bSort);
603
604 ee_sum_nsteps = 0;
605 ee_sum_nsum = 0;
606 snew(ee_sum, nremax)(ee_sum) = save_calloc("ee_sum", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 606, (nremax), sizeof(*(ee_sum)));
607
608 snew(fr, 1)(fr) = save_calloc("fr", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 608, (1), sizeof(*(fr)));
609 snew(fro, 1)(fro) = save_calloc("fro", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 609, (1), sizeof(*(fro)));
610 fro->t = -1e20;
611 fro->nre = nre;
612 snew(fro->ener, nremax)(fro->ener) = save_calloc("fro->ener", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 612, (nremax), sizeof(*(fro->ener)));
613
614 noutfr = 0;
615 bFirst = TRUE1;
616
617 last_t = fro->t;
618 for (f = 0; f < nfile; f++)
619 {
620 bNewFile = TRUE1;
621 bNewOutput = TRUE1;
622 in = open_enx(fnms[f], "r");
623 enm = NULL((void*)0);
624 do_enxnms(in, &this_nre, &enm);
625 if (f == 0)
626 {
627 if (scalefac != 1)
628 {
629 set = select_it(nre, enm, &nset);
630 }
631
632 /* write names to the output file */
633 out = open_enx(opt2fn("-o", NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "w");
634 do_enxnms(out, &nre, &enm);
635 }
636
637 /* start reading from the next file */
638 while ((fro->t <= (settime[f+1] + GMX_REAL_EPS5.96046448E-08)) &&
639 do_enx(in, fr))
640 {
641 if (bNewFile)
642 {
643 startstep_file = fr->step;
644 tadjust = settime[f] - fr->t;
645 if (cont_type[f+1] == TIME_LAST2)
646 {
647 settime[f+1] = readtime[f+1]-readtime[f]+settime[f];
648 cont_type[f+1] = TIME_EXPLICIT0;
649 }
650 bNewFile = FALSE0;
651 }
652
653 if (tadjust + fr->t <= last_t)
654 {
655 /* Skip this frame, since we already have it / past it */
656 if (debug)
657 {
658 fprintf(debug, "fr->step %s, fr->t %.4f\n",
659 gmx_step_str(fr->step, buf), fr->t);
660 fprintf(debug, "tadjust %12.6e + fr->t %12.6e <= t %12.6e\n",
661 tadjust, fr->t, last_t);
662 }
663 continue;
664 }
665
666 fro->step = lastfilestep + fr->step - startstep_file;
667 fro->t = tadjust + fr->t;
668
669 bWrite = ((begin < 0 || (begin >= 0 && (fro->t >= begin-GMX_REAL_EPS5.96046448E-08))) &&
670 (end < 0 || (end >= 0 && (fro->t <= end +GMX_REAL_EPS5.96046448E-08))) &&
671 (fro->t <= settime[f+1]+0.5*timestep));
672
673 if (debug)
674 {
675 fprintf(debug,
676 "fr->step %s, fr->t %.4f, fro->step %s fro->t %.4f, w %d\n",
677 gmx_step_str(fr->step, buf), fr->t,
678 gmx_step_str(fro->step, buf2), fro->t, bWrite);
679 }
680
681 if (bError)
682 {
683 if ((end > 0) && (fro->t > end+GMX_REAL_EPS5.96046448E-08))
684 {
685 f = nfile;
686 break;
687 }
688 }
689
690 if (fro->t >= begin-GMX_REAL_EPS5.96046448E-08)
691 {
692 if (bFirst)
693 {
694 bFirst = FALSE0;
695 startstep = fr->step;
696 }
697 if (bWrite)
698 {
699 update_ee_sum(nre, &ee_sum_step, &ee_sum_nsteps, &ee_sum_nsum, ee_sum,
700 fr, fro->step);
701 }
702 }
703
704 /* determine if we should write it */
705 if (bWrite && (delta_t == 0 || bRmod(fro->t, toffset, delta_t)bRmod_fd(fro->t, toffset, delta_t, 0)))
706 {
707 laststep = fro->step;
708 last_t = fro->t;
709 if (bNewOutput)
710 {
711 bNewOutput = FALSE0;
712 fprintf(stderrstderr, "\nContinue writing frames from t=%g, step=%s\n",
713 fro->t, gmx_step_str(fro->step, buf));
714 }
715
716 /* Copy the energies */
717 for (i = 0; i < nre; i++)
718 {
719 fro->ener[i].e = fr->ener[i].e;
720 }
721
722 fro->nsteps = ee_sum_nsteps;
723 fro->dt = fr->dt;
724
725 if (ee_sum_nsum <= 1)
726 {
727 fro->nsum = 0;
728 }
729 else
730 {
731 fro->nsum = gmx_int64_to_int(ee_sum_nsum,
732 "energy average summation");
733 /* Copy the energy sums */
734 for (i = 0; i < nre; i++)
735 {
736 fro->ener[i].esum = ee_sum[i].esum;
737 fro->ener[i].eav = ee_sum[i].eav;
738 }
739 }
740 /* We wrote the energies, so reset the counts */
741 ee_sum_nsteps = 0;
742 ee_sum_nsum = 0;
743
744 if (scalefac != 1)
745 {
746 for (kkk = 0; kkk < nset; kkk++)
747 {
748 fro->ener[set[kkk]].e *= scalefac;
749 if (fro->nsum > 0)
750 {
751 fro->ener[set[kkk]].eav *= scalefac*scalefac;
752 fro->ener[set[kkk]].esum *= scalefac;
753 }
754 }
755 }
756 /* Copy restraint stuff */
757 /*fro->ndisre = fr->ndisre;
758 fro->disre_rm3tav = fr->disre_rm3tav;
759 fro->disre_rt = fr->disre_rt;*/
760 fro->nblock = fr->nblock;
761 /*fro->nr = fr->nr;*/
762 fro->block = fr->block;
763
764 /* check if we have blocks with delta_h data and are throwing
765 away data */
766 if (fro->nblock > 0)
767 {
768 if (remove_dh)
769 {
770 int i;
771 if (!blocks || nblocks_alloc < fr->nblock)
772 {
773 /* we pre-allocate the blocks */
774 nblocks_alloc = fr->nblock;
775 snew(blocks, nblocks_alloc)(blocks) = save_calloc("blocks", "/home/alexxy/Develop/gromacs/src/gromacs/gmxana/gmx_eneconv.c"
, 775, (nblocks_alloc), sizeof(*(blocks)));
776 }
777 nblocks = 0; /* number of blocks so far */
778
779 for (i = 0; i < fr->nblock; i++)
780 {
781 if ( (fr->block[i].id != enxDHCOLL) &&
782 (fr->block[i].id != enxDH) &&
783 (fr->block[i].id != enxDHHIST) )
784 {
785 /* copy everything verbatim */
786 blocks[nblocks] = fr->block[i];
787 nblocks++;
788 }
789 }
790 /* now set the block pointer to the new blocks */
791 fro->nblock = nblocks;
792 fro->block = blocks;
793 }
794 else if (delta_t > 0)
795 {
796 if (!warned_about_dh)
797 {
798 for (i = 0; i < fr->nblock; i++)
799 {
800 if (fr->block[i].id == enxDH ||
801 fr->block[i].id == enxDHHIST)
802 {
803 int size;
804 if (fr->block[i].id == enxDH)
805 {
806 size = fr->block[i].sub[2].nr;
807 }
808 else
809 {
810 size = fr->nsteps;
811 }
812 if (size > 0)
813 {
814 printf("\nWARNING: %s contains delta H blocks or histograms for which\n"
815 " some data is thrown away on a block-by-block basis, where each block\n"
816 " contains up to %d samples.\n"
817 " This is almost certainly not what you want.\n"
818 " Use the -rmdh option to throw all delta H samples away.\n"
819 " Use g_energy -odh option to extract these samples.\n",
820 fnms[f], size);
821 warned_about_dh = TRUE1;
822 break;
823 }
824 }
825 }
826 }
827 }
828 }
829
830 do_enx(out, fro);
831 if (noutfr % 1000 == 0)
832 {
833 fprintf(stderrstderr, "Writing frame time %g ", fro->t);
834 }
835 noutfr++;
836 }
837 }
838 if (f == nfile)
839 {
840 f--;
841 }
842 printf("\nLast step written from %s: t %g, step %s\n",
843 fnms[f], last_t, gmx_step_str(laststep, buf));
844 lastfilestep = laststep;
845
846 /* set the next time from the last in previous file */
847 if (cont_type[f+1] == TIME_CONTINUE1)
848 {
849 settime[f+1] = fro->t;
850 /* in this case we have already written the last frame of
851 * previous file, so update begin to avoid doubling it
852 * with the start of the next file
853 */
854 begin = fro->t+0.5*timestep;
855 /* cont_type[f+1]==TIME_EXPLICIT; */
856 }
857
858 if ((fro->t < end) && (f < nfile-1) &&
859 (fro->t < settime[f+1]-1.5*timestep))
860 {
861 fprintf(stderrstderr,
862 "\nWARNING: There might be a gap around t=%g\n", fro->t);
863 }
864
865 /* move energies to lastee */
866 close_enx(in);
867 free_enxnms(this_nre, enm);
868
869 fprintf(stderrstderr, "\n");
870 }
871 if (noutfr == 0)
872 {
873 fprintf(stderrstderr, "No frames written.\n");
874 }
875 else
876 {
877 fprintf(stderrstderr, "Last frame written was at step %s, time %f\n",
878 gmx_step_str(fro->step, buf), fro->t);
879 fprintf(stderrstderr, "Wrote %d frames\n", noutfr);
880 }
881
882 return 0;
883}