?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge branch 'master' of git://git.gromacs.org/gromacs
[gromacs/adressmacs.git] / src / tools / gmx_eneconv.c
blob09071ff85ad9024edf039ed7d7bf27cf229ce0bb
1 /*
2 *
3 * This source code is part of
4 *
5 * G R O M A C S
6 *
7 * GROningen MAchine for Chemical Simulations
8 *
9 * VERSION 3.2.0
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
14
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
19 *
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
26 *
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
29 *
30 * For more info, check our website at http://www.gromacs.org
31 *
32 * And Hey:
33 * Green Red Orange Magenta Azure Cyan Skyblue
34 */
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38 #include <string.h>
39 #include <math.h>
40
41 #include "string2.h"
42 #include "typedefs.h"
43 #include "smalloc.h"
44 #include "statutil.h"
45 #include "disre.h"
46 #include "names.h"
47 #include "copyrite.h"
48 #include "macros.h"
49 #include "gmx_fatal.h"
50 #include "enxio.h"
51 #include "vec.h"
52 #include "gmx_ana.h"
53
54 #define TIME_EXPLICIT 0
55 #define TIME_CONTINUE 1
56 #define TIME_LAST 2
57 #ifndef FLT_MAX
58 #define FLT_MAX 1e36
59 #endif
60
61 static int *select_it(int nre,gmx_enxnm_t *nm,int *nset)
62 {
63 gmx_bool *bE;
64 int n,k,j,i;
65 int *set;
66 gmx_bool bVerbose = TRUE;
67
68 if ((getenv("VERBOSE")) != NULL)
69 bVerbose = FALSE;
70
71 fprintf(stderr,"Select the terms you want to scale from the following list\n");
72 fprintf(stderr,"End your selection with 0\n");
73
74 if ( bVerbose ) {
75 for(k=0; (k<nre); ) {
76 for(j=0; (j<4) && (k<nre); j++,k++)
77 fprintf(stderr," %3d=%14s",k+1,nm[k].name);
78 fprintf(stderr,"\n");
79 }
80 }
81
82 snew(bE,nre);
83 do {
84 if(1 != scanf("%d",&n))
85 {
86 gmx_fatal(FARGS,"Cannot read energy term");
87 }
88 if ((n>0) && (n<=nre))
89 bE[n-1]=TRUE;
90 } while (n != 0);
91
92 snew(set,nre);
93 for(i=(*nset)=0; (i<nre); i++)
94 if (bE[i])
95 set[(*nset)++]=i;
96
97 sfree(bE);
98
99 return set;
100 }
101
102 static gmx_bool same_time(real t1,real t2)
103 {
104 const real tol=1e-5;
105
106 return (fabs(t1-t2) < tol);
107 }
108
109
110 gmx_bool bRgt(double a,double b)
111 {
112 double tol = 1e-6;
113
114 if ( a > (b - tol*(a+b)) )
115 return TRUE;
116 else
117 return FALSE;
118 }
119
120 static void sort_files(char **fnms,real *settime,int nfile)
121 {
122 int i,j,minidx;
123 real timeswap;
124 char *chptr;
125
126 for(i=0;i<nfile;i++) {
127 minidx=i;
128 for(j=i+1;j<nfile;j++) {
129 if(settime[j]<settime[minidx])
130 minidx=j;
131 }
132 if(minidx!=i) {
133 timeswap=settime[i];
134 settime[i]=settime[minidx];
135 settime[minidx]=timeswap;
136 chptr=fnms[i];
137 fnms[i]=fnms[minidx];
138 fnms[minidx]=chptr;
139 }
140 }
141 }
142
143
144 static int scan_ene_files(char **fnms, int nfiles,
145 real *readtime, real *timestep, int *nremax)
146 {
147 /* Check number of energy terms and start time of all files */
148 int f,i,nre,nremin=0,nresav=0;
149 ener_file_t in;
150 real t1,t2;
151 char inputstring[STRLEN];
152 gmx_enxnm_t *enm;
153 t_enxframe *fr;
154
155 snew(fr,1);
156
157 for(f=0; f<nfiles; f++) {
158 in = open_enx(fnms[f],"r");
159 enm = NULL;
160 do_enxnms(in,&nre,&enm);
161
162 if (f == 0) {
163 nresav = nre;
164 nremin = nre;
165 *nremax = nre;
166 do_enx(in,fr);
167 t1 = fr->t;
168 do_enx(in,fr);
169 t2 = fr->t;
170 *timestep=t2-t1;
171 readtime[f]=t1;
172 close_enx(in);
173 } else {
174 nremin = min(nremin,fr->nre);
175 *nremax = max(*nremax,fr->nre);
176 if (nre != nresav) {
177 fprintf(stderr,
178 "Energy files don't match, different number of energies:\n"
179 " %s: %d\n %s: %d\n",fnms[f-1],nresav,fnms[f],fr->nre);
180 fprintf(stderr,
181 "\nContinue conversion using only the first %d terms (n/y)?\n"
182 "(you should be sure that the energy terms match)\n",nremin);
183 if(NULL==fgets(inputstring,STRLEN-1,stdin))
184 {
185 gmx_fatal(FARGS,"Error reading user input");
186 }
187 if (inputstring[0]!='y' && inputstring[0]!='Y') {
188 fprintf(stderr,"Will not convert\n");
189 exit(0);
190 }
191 nresav = fr->nre;
192 }
193 do_enx(in,fr);
194 readtime[f] = fr->t;
195 close_enx(in);
196 }
197 fprintf(stderr,"\n");
198 free_enxnms(nre,enm);
199 }
200
201 free_enxframe(fr);
202 sfree(fr);
203
204 return nremin;
205 }
206
207
208 static void edit_files(char **fnms,int nfiles,real *readtime,
209 real *settime,int *cont_type,gmx_bool bSetTime,gmx_bool bSort)
210 {
211 int i;
212 gmx_bool ok;
213 char inputstring[STRLEN],*chptr;
214
215 if(bSetTime) {
216 if(nfiles==1)
217 fprintf(stderr,"\n\nEnter the new start time:\n\n");
218 else
219 fprintf(stderr,"\n\nEnter the new start time for each file.\n"
220 "There are two special options, both disables sorting:\n\n"
221 "c (continue) - The start time is taken from the end\n"
222 "of the previous file. Use it when your continuation run\n"
223 "restarts with t=0 and there is no overlap.\n\n"
224 "l (last) - The time in this file will be changed the\n"
225 "same amount as in the previous. Use it when the time in the\n"
226 "new run continues from the end of the previous one,\n"
227 "since this takes possible overlap into account.\n\n");
228
229 fprintf(stderr," File Current start New start\n"
230 "---------------------------------------------------------\n");
231
232 for(i=0;i<nfiles;i++) {
233 fprintf(stderr,"%25s %10.3f ",fnms[i],readtime[i]);
234 ok=FALSE;
235 do {
236 if(NULL==fgets(inputstring,STRLEN-1,stdin))
237 {
238 gmx_fatal(FARGS,"Error reading user input");
239 }
240 inputstring[strlen(inputstring)-1]=0;
241
242 if(inputstring[0]=='c' || inputstring[0]=='C') {
243 cont_type[i]=TIME_CONTINUE;
244 bSort=FALSE;
245 ok=TRUE;
246 settime[i]=FLT_MAX;
247 }
248 else if(inputstring[0]=='l' ||
249 inputstring[0]=='L') {
250 cont_type[i]=TIME_LAST;
251 bSort=FALSE;
252 ok=TRUE;
253 settime[i]=FLT_MAX;
254 }
255 else {
256 settime[i]=strtod(inputstring,&chptr);
257 if(chptr==inputstring) {
258 fprintf(stderr,"Try that again: ");
259 }
260 else {
261 cont_type[i]=TIME_EXPLICIT;
262 ok=TRUE;
263 }
264 }
265 } while (!ok);
266 }
267 if(cont_type[0]!=TIME_EXPLICIT) {
268 cont_type[0]=TIME_EXPLICIT;
269 settime[0]=0;
270 }
271 }
272 else
273 for(i=0;i<nfiles;i++)
274 settime[i]=readtime[i];
275
276 if(bSort && (nfiles>1))
277 sort_files(fnms,settime,nfiles);
278 else
279 fprintf(stderr,"Sorting disabled.\n");
280
281
282 /* Write out the new order and start times */
283 fprintf(stderr,"\nSummary of files and start times used:\n\n"
284 " File Start time\n"
285 "-----------------------------------------\n");
286 for(i=0;i<nfiles;i++)
287 switch(cont_type[i]) {
288 case TIME_EXPLICIT:
289 fprintf(stderr,"%25s %10.3f\n",fnms[i],settime[i]);
290 break;
291 case TIME_CONTINUE:
292 fprintf(stderr,"%25s Continue from end of last file\n",fnms[i]);
293 break;
294 case TIME_LAST:
295 fprintf(stderr,"%25s Change by same amount as last file\n",fnms[i]);
296 break;
297 }
298 fprintf(stderr,"\n");
299
300 settime[nfiles]=FLT_MAX;
301 cont_type[nfiles]=TIME_EXPLICIT;
302 readtime[nfiles]=FLT_MAX;
303 }
304
305
306 static void copy_ee(t_energy *src, t_energy *dst, int nre)
307 {
308 int i;
309
310 for(i=0;i<nre;i++) {
311 dst[i].e=src[i].e;
312 dst[i].esum=src[i].esum;
313 dst[i].eav=src[i].eav;
314 }
315 }
316
317
318 static void remove_last_eeframe(t_energy *lastee, gmx_large_int_t laststep,
319 t_energy *ee, int nre)
320 {
321 int i;
322 gmx_large_int_t p=laststep+1;
323 double sigmacorr;
324
325 for(i=0;i<nre;i++) {
326 lastee[i].esum-=ee[i].e;
327 sigmacorr=lastee[i].esum-(p-1)*ee[i].e;
328 lastee[i].eav-=(sigmacorr*sigmacorr)/((p-1)*p);
329 }
330 }
331
332
333
334 static void update_ee(t_energy *lastee,gmx_large_int_t laststep,
335 t_energy *startee,gmx_large_int_t startstep,
336 t_energy *ee, int step,
337 t_energy *outee, int nre)
338 {
339 int i;
340 double sigmacorr,nom,denom;
341 double prestart_esum;
342 double prestart_sigma;
343
344 for(i=0;i<nre;i++) {
345 outee[i].e=ee[i].e;
346 /* add statistics from earlier file if present */
347 if(laststep>0) {
348 outee[i].esum=lastee[i].esum+ee[i].esum;
349 nom=(lastee[i].esum*(step+1)-ee[i].esum*(laststep));
350 denom=((step+1.0)*(laststep)*(step+1.0+laststep));
351 sigmacorr=nom*nom/denom;
352 outee[i].eav=lastee[i].eav+ee[i].eav+sigmacorr;
353 }
354 else {
355 /* otherwise just copy to output */
356 outee[i].esum=ee[i].esum;
357 outee[i].eav=ee[i].eav;
358 }
359
360 /* if we didnt start to write at the first frame
361 * we must compensate the statistics for this
362 * there are laststep frames in the earlier file
363 * and step+1 frames in this one.
364 */
365 if(startstep>0) {
366 gmx_large_int_t q=laststep+step;
367 gmx_large_int_t p=startstep+1;
368 prestart_esum=startee[i].esum-startee[i].e;
369 sigmacorr=prestart_esum-(p-1)*startee[i].e;
370 prestart_sigma=startee[i].eav-
371 sigmacorr*sigmacorr/(p*(p-1));
372 sigmacorr=prestart_esum/(p-1)-
373 outee[i].esum/(q);
374 outee[i].esum-=prestart_esum;
375 if (q-p+1 > 0)
376 outee[i].eav=outee[i].eav-prestart_sigma-
377 sigmacorr*sigmacorr*((p-1)*q)/(q-p+1);
378 }
379
380 if((outee[i].eav/(laststep+step+1))<(GMX_REAL_EPS))
381 outee[i].eav=0;
382 }
383 }
384
385
386 static void update_last_ee(t_energy *lastee, gmx_large_int_t laststep,
387 t_energy *ee,gmx_large_int_t step,int nre)
388 {
389 t_energy *tmp;
390 snew(tmp,nre);
391 update_ee(lastee,laststep,NULL,0,ee,step,tmp,nre);
392 copy_ee(tmp,lastee,nre);
393 sfree(tmp);
394 }
395
396 static void update_ee_sum(int nre,
397 gmx_large_int_t *ee_sum_step,
398 gmx_large_int_t *ee_sum_nsteps,
399 gmx_large_int_t *ee_sum_nsum,
400 t_energy *ee_sum,
401 t_enxframe *fr,int out_step)
402 {
403 gmx_large_int_t nsteps,nsum,fr_nsum;
404 int i;
405
406 nsteps = *ee_sum_nsteps;
407 nsum = *ee_sum_nsum;
408
409 fr_nsum = fr->nsum;
410 if (fr_nsum == 0) {
411 fr_nsum = 1;
412 }
413
414 if (nsteps == 0) {
415 if (fr_nsum == 1) {
416 for(i=0;i<nre;i++) {
417 ee_sum[i].esum = fr->ener[i].e;
418 ee_sum[i].eav = 0;
419 }
420 } else {
421 for(i=0;i<nre;i++) {
422 ee_sum[i].esum = fr->ener[i].esum;
423 ee_sum[i].eav = fr->ener[i].eav;
424 }
425 }
426 nsteps = fr->nsteps;
427 nsum = fr_nsum;
428 } else if (out_step + *ee_sum_nsum - *ee_sum_step == nsteps + fr->nsteps) {
429 if (fr_nsum == 1) {
430 for(i=0;i<nre;i++) {
431 ee_sum[i].eav +=
432 dsqr(ee_sum[i].esum/nsum
433 - (ee_sum[i].esum + fr->ener[i].e)/(nsum + 1))*nsum*(nsum + 1);
434 ee_sum[i].esum += fr->ener[i].e;
435 }
436 } else {
437 for(i=0; i<fr->nre; i++) {
438 ee_sum[i].eav +=
439 fr->ener[i].eav +
440 dsqr(ee_sum[i].esum/nsum
441 - (ee_sum[i].esum + fr->ener[i].esum)/(nsum + fr->nsum))*
442 nsum*(nsum + fr->nsum)/(double)fr->nsum;
443 ee_sum[i].esum += fr->ener[i].esum;
444 }
445 }
446 nsteps += fr->nsteps;
447 nsum += fr_nsum;
448 } else {
449 if (fr->nsum != 0) {
450 fprintf(stderr,"\nWARNING: missing energy sums at time %f\n",fr->t);
451 }
452 nsteps = 0;
453 nsum = 0;
454 }
455
456 *ee_sum_step = out_step;
457 *ee_sum_nsteps = nsteps;
458 *ee_sum_nsum = nsum;
459 }
460
461 int gmx_eneconv(int argc,char *argv[])
462 {
463 const char *desc[] = {
464 "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[BR]",
465 "Concatenates several energy files in sorted order.",
466 "In case of double time frames the one",
467 "in the later file is used. By specifying [TT]-settime[tt] you will be",
468 "asked for the start time of each file. The input files are taken",
469 "from the command line,",
470 "such that the command [TT]eneconv -o fixed.edr *.edr[tt] should do",
471 "the trick. [PAR]",
472 "With [IT]one file[it] specified for [TT]-f[tt]:[BR]",
473 "Reads one energy file and writes another, applying the [TT]-dt[tt],",
474 "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
475 "converting to a different format if necessary (indicated by file",
476 "extentions).[PAR]",
477 "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
478 "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
479 };
480 const char *bugs[] = {
481 "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."
482 };
483 ener_file_t in=NULL, out=NULL;
484 gmx_enxnm_t *enm=NULL;
485 #if 0
486 ener_file_t in,out=NULL;
487 gmx_enxnm_t *enm=NULL;
488 #endif
489 t_enxframe *fr,*fro;
490 gmx_large_int_t ee_sum_step=0,ee_sum_nsteps,ee_sum_nsum;
491 t_energy *ee_sum;
492 gmx_large_int_t lastfilestep,laststep,startstep,startstep_file=0;
493 int noutfr;
494 int nre,nremax,this_nre,nfile,f,i,j,kkk,nset,*set=NULL;
495 double last_t;
496 char **fnms;
497 real *readtime,*settime,timestep,t1,tadjust;
498 char inputstring[STRLEN],*chptr,buf[22],buf2[22],buf3[22];
499 gmx_bool ok;
500 int *cont_type;
501 gmx_bool bNewFile,bFirst,bNewOutput;
502 output_env_t oenv;
503 gmx_bool warned_about_dh=FALSE;
504 t_enxblock *blocks=NULL;
505 int nblocks=0;
506 int nblocks_alloc=0;
507
508 t_filenm fnm[] = {
509 { efEDR, "-f", NULL, ffRDMULT },
510 { efEDR, "-o", "fixed", ffWRITE },
511 };
512
513 #define NFILE asize(fnm)
514 gmx_bool bWrite;
515 static real delta_t=0.0, toffset=0,scalefac=1;
516 static gmx_bool bSetTime=FALSE;
517 static gmx_bool bSort=TRUE,bError=TRUE;
518 static real begin=-1;
519 static real end=-1;
520 gmx_bool remove_dh=FALSE;
521
522 t_pargs pa[] = {
523 { "-b", FALSE, etREAL, {&begin},
524 "First time to use"},
525 { "-e", FALSE, etREAL, {&end},
526 "Last time to use"},
527 { "-dt", FALSE, etREAL, {&delta_t},
528 "Only write out frame when t MOD dt = offset" },
529 { "-offset", FALSE, etREAL, {&toffset},
530 "Time offset for -dt option" },
531 { "-settime", FALSE, etBOOL, {&bSetTime},
532 "Change starting time interactively" },
533 { "-sort", FALSE, etBOOL, {&bSort},
534 "Sort energy files (not frames)"},
535 { "-rmdh", FALSE, etBOOL, {&remove_dh},
536 "Remove free energy block data" },
537 { "-scalefac", FALSE, etREAL, {&scalefac},
538 "Multiply energy component by this factor" },
539 { "-error", FALSE, etBOOL, {&bError},
540 "Stop on errors in the file" }
541 };
542
543 CopyRight(stderr,argv[0]);
544 parse_common_args(&argc,argv,PCA_BE_NICE,NFILE,fnm,asize(pa),
545 pa,asize(desc),desc,asize(bugs),bugs,&oenv);
546 tadjust = 0;
547 nremax = 0;
548 nset = 0;
549 timestep = 0.0;
550 snew(fnms,argc);
551 nfile=0;
552 lastfilestep=0;
553 laststep=startstep=0;
554
555 nfile = opt2fns(&fnms,"-f",NFILE,fnm);
556
557 if (!nfile)
558 gmx_fatal(FARGS,"No input files!");
559
560 snew(settime,nfile+1);
561 snew(readtime,nfile+1);
562 snew(cont_type,nfile+1);
563
564 nre=scan_ene_files(fnms,nfile,readtime,&timestep,&nremax);
565 edit_files(fnms,nfile,readtime,settime,cont_type,bSetTime,bSort);
566
567 ee_sum_nsteps = 0;
568 ee_sum_nsum = 0;
569 snew(ee_sum,nremax);
570
571 snew(fr,1);
572 snew(fro,1);
573 fro->t = -1e20;
574 fro->nre = nre;
575 snew(fro->ener,nremax);
576
577 noutfr=0;
578 bFirst=TRUE;
579
580 last_t = fro->t;
581 for(f=0;f<nfile;f++) {
582 bNewFile=TRUE;
583 bNewOutput=TRUE;
584 in=open_enx(fnms[f],"r");
585 enm = NULL;
586 do_enxnms(in,&this_nre,&enm);
587 if (f == 0) {
588 if (scalefac != 1)
589 set = select_it(nre,enm,&nset);
590
591 /* write names to the output file */
592 out=open_enx(opt2fn("-o",NFILE,fnm),"w");
593 do_enxnms(out,&nre,&enm);
594 }
595
596 /* start reading from the next file */
597 while ((fro->t <= (settime[f+1] + GMX_REAL_EPS)) &&
598 do_enx(in,fr)) {
599 if(bNewFile) {
600 startstep_file = fr->step;
601 tadjust = settime[f] - fr->t;
602 if(cont_type[f+1]==TIME_LAST) {
603 settime[f+1] = readtime[f+1]-readtime[f]+settime[f];
604 cont_type[f+1] = TIME_EXPLICIT;
605 }
606 bNewFile = FALSE;
607 }
608
609 if (tadjust + fr->t <= last_t) {
610 /* Skip this frame, since we already have it / past it */
611 if (debug) {
612 fprintf(debug,"fr->step %s, fr->t %.4f\n",
613 gmx_step_str(fr->step,buf),fr->t);
614 fprintf(debug,"tadjust %12.6e + fr->t %12.6e <= t %12.6e\n",
615 tadjust,fr->t,last_t);
616 }
617 continue;
618 }
619
620 fro->step = lastfilestep + fr->step - startstep_file;
621 fro->t = tadjust + fr->t;
622
623 bWrite = ((begin<0 || (begin>=0 && (fro->t >= begin-GMX_REAL_EPS))) &&
624 (end <0 || (end >=0 && (fro->t <= end +GMX_REAL_EPS))) &&
625 (fro->t <= settime[f+1]+0.5*timestep));
626
627 if (debug) {
628 fprintf(debug,
629 "fr->step %s, fr->t %.4f, fro->step %s fro->t %.4f, w %d\n",
630 gmx_step_str(fr->step,buf),fr->t,
631 gmx_step_str(fro->step,buf2),fro->t,bWrite);
632 }
633
634 if (bError)
635 if ((end > 0) && (fro->t > end+GMX_REAL_EPS)) {
636 f = nfile;
637 break;
638 }
639
640 if (fro->t >= begin-GMX_REAL_EPS) {
641 if (bFirst) {
642 bFirst = FALSE;
643 startstep = fr->step;
644 }
645 if (bWrite) {
646 update_ee_sum(nre,&ee_sum_step,&ee_sum_nsteps,&ee_sum_nsum,ee_sum,
647 fr,fro->step);
648 }
649 }
650
651 /* determine if we should write it */
652 if (bWrite && (delta_t==0 || bRmod(fro->t,toffset,delta_t))) {
653 laststep = fro->step;
654 last_t = fro->t;
655 if(bNewOutput) {
656 bNewOutput=FALSE;
657 fprintf(stderr,"\nContinue writing frames from t=%g, step=%s\n",
658 fro->t,gmx_step_str(fro->step,buf));
659 }
660
661 /* Copy the energies */
662 for(i=0; i<nre; i++) {
663 fro->ener[i].e = fr->ener[i].e;
664 }
665
666 fro->nsteps = ee_sum_nsteps;
667 fro->dt = fr->dt;
668
669 if (ee_sum_nsum <= 1) {
670 fro->nsum = 0;
671 } else {
672 fro->nsum = gmx_large_int_to_int(ee_sum_nsum,
673 "energy average summation");
674 /* Copy the energy sums */
675 for(i=0; i<nre; i++) {
676 fro->ener[i].esum = ee_sum[i].esum;
677 fro->ener[i].eav = ee_sum[i].eav;
678 }
679 }
680 /* We wrote the energies, so reset the counts */
681 ee_sum_nsteps = 0;
682 ee_sum_nsum = 0;
683
684 if (scalefac != 1) {
685 for(kkk=0; kkk<nset; kkk++) {
686 fro->ener[set[kkk]].e *= scalefac;
687 if (fro->nsum > 0) {
688 fro->ener[set[kkk]].eav *= scalefac*scalefac;
689 fro->ener[set[kkk]].esum *= scalefac;
690 }
691 }
692 }
693 /* Copy restraint stuff */
694 /*fro->ndisre = fr->ndisre;
695 fro->disre_rm3tav = fr->disre_rm3tav;
696 fro->disre_rt = fr->disre_rt;*/
697 fro->nblock = fr->nblock;
698 /*fro->nr = fr->nr;*/
699 fro->block = fr->block;
700
701 /* check if we have blocks with delta_h data and are throwing
702 away data */
703 if (fro->nblock > 0)
704 {
705 if (remove_dh)
706 {
707 int i;
708 if (!blocks || nblocks_alloc < fr->nblock)
709 {
710 /* we pre-allocate the blocks */
711 nblocks_alloc=fr->nblock;
712 snew(blocks, nblocks_alloc);
713 }
714 nblocks=0; /* number of blocks so far */
715
716 for(i=0;i<fr->nblock;i++)
717 {
718 if ( (fr->block[i].id != enxDHCOLL) &&
719 (fr->block[i].id != enxDH) &&
720 (fr->block[i].id != enxDHHIST) )
721 {
722 /* copy everything verbatim */
723 blocks[nblocks] = fr->block[i];
724 nblocks++;
725 }
726 }
727 /* now set the block pointer to the new blocks */
728 fro->nblock = nblocks;
729 fro->block = blocks;
730 }
731 else if (delta_t > 0)
732 {
733 if (!warned_about_dh)
734 {
735 for(i=0;i<fr->nblock;i++)
736 {
737 if (fr->block[i].id == enxDH ||
738 fr->block[i].id == enxDHHIST)
739 {
740 int size;
741 if (fr->block[i].id == enxDH )
742 {
743 size=fr->block[i].sub[2].nr;
744 }
745 else
746 {
747 size=fr->nsteps;
748 }
749 if (size>0)
750 {
751 printf("\nWARNING: %s contains delta H blocks or histograms for which\n"
752 " some data is thrown away on a block-by-block basis, where each block\n"
753 " contains up to %d samples.\n"
754 " This is almost certainly not what you want.\n"
755 " Use the -rmdh option to throw all delta H samples away.\n"
756 " Use g_energy -odh option to extract these samples.\n",
757 fnms[f], size);
758 warned_about_dh=TRUE;
759 break;
760 }
761 }
762 }
763 }
764 }
765 }
766
767 do_enx(out,fro);
768 if (noutfr % 1000 == 0)
769 fprintf(stderr,"Writing frame time %g ",fro->t);
770 noutfr++;
771 }
772 }
773
774 printf("\nLast step written from %s: t %g, step %s\n",
775 fnms[f],last_t,gmx_step_str(laststep,buf));
776 lastfilestep = laststep;
777
778 /* set the next time from the last in previous file */
779 if (cont_type[f+1]==TIME_CONTINUE) {
780 settime[f+1] = fro->t;
781 /* in this case we have already written the last frame of
782 * previous file, so update begin to avoid doubling it
783 * with the start of the next file
784 */
785 begin = fro->t+0.5*timestep;
786 /* cont_type[f+1]==TIME_EXPLICIT; */
787 }
788
789 if ((fro->t < end) && (f < nfile-1) &&
790 (fro->t < settime[f+1]-1.5*timestep))
791 fprintf(stderr,
792 "\nWARNING: There might be a gap around t=%g\n",fro->t);
793
794 /* move energies to lastee */
795 close_enx(in);
796 free_enxnms(this_nre,enm);
797
798 fprintf(stderr,"\n");
799 }
800 if (noutfr == 0)
801 fprintf(stderr,"No frames written.\n");
802 else {
803 fprintf(stderr,"Last frame written was at step %s, time %f\n",
804 gmx_step_str(fro->step,buf),fro->t);
805 fprintf(stderr,"Wrote %d frames\n",noutfr);
806 }
807
808 thanx(stderr);
809 return 0;
810 }
got merge a long time ago