3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
49 #include "gmx_fatal.h"
54 #define TIME_EXPLICIT 0
55 #define TIME_CONTINUE 1
61 static int *select_it(int nre
,gmx_enxnm_t
*nm
,int *nset
)
66 gmx_bool bVerbose
= TRUE
;
68 if ((getenv("VERBOSE")) != NULL
)
71 fprintf(stderr
,"Select the terms you want to scale from the following list\n");
72 fprintf(stderr
,"End your selection with 0\n");
76 for(j
=0; (j
<4) && (k
<nre
); j
++,k
++)
77 fprintf(stderr
," %3d=%14s",k
+1,nm
[k
].name
);
84 if(1 != scanf("%d",&n
))
86 gmx_fatal(FARGS
,"Cannot read energy term");
88 if ((n
>0) && (n
<=nre
))
93 for(i
=(*nset
)=0; (i
<nre
); i
++)
102 static void sort_files(char **fnms
,real
*settime
,int nfile
)
108 for(i
=0;i
<nfile
;i
++) {
110 for(j
=i
+1;j
<nfile
;j
++) {
111 if(settime
[j
]<settime
[minidx
])
116 settime
[i
]=settime
[minidx
];
117 settime
[minidx
]=timeswap
;
119 fnms
[i
]=fnms
[minidx
];
126 static int scan_ene_files(char **fnms
, int nfiles
,
127 real
*readtime
, real
*timestep
, int *nremax
)
129 /* Check number of energy terms and start time of all files */
130 int f
,i
,nre
,nremin
=0,nresav
=0;
133 char inputstring
[STRLEN
];
139 for(f
=0; f
<nfiles
; f
++) {
140 in
= open_enx(fnms
[f
],"r");
142 do_enxnms(in
,&nre
,&enm
);
156 nremin
= min(nremin
,fr
->nre
);
157 *nremax
= max(*nremax
,fr
->nre
);
160 "Energy files don't match, different number of energies:\n"
161 " %s: %d\n %s: %d\n",fnms
[f
-1],nresav
,fnms
[f
],fr
->nre
);
163 "\nContinue conversion using only the first %d terms (n/y)?\n"
164 "(you should be sure that the energy terms match)\n",nremin
);
165 if(NULL
==fgets(inputstring
,STRLEN
-1,stdin
))
167 gmx_fatal(FARGS
,"Error reading user input");
169 if (inputstring
[0]!='y' && inputstring
[0]!='Y') {
170 fprintf(stderr
,"Will not convert\n");
179 fprintf(stderr
,"\n");
180 free_enxnms(nre
,enm
);
190 static void edit_files(char **fnms
,int nfiles
,real
*readtime
,
191 real
*settime
,int *cont_type
,gmx_bool bSetTime
,gmx_bool bSort
)
195 char inputstring
[STRLEN
],*chptr
;
199 fprintf(stderr
,"\n\nEnter the new start time:\n\n");
201 fprintf(stderr
,"\n\nEnter the new start time for each file.\n"
202 "There are two special options, both disables sorting:\n\n"
203 "c (continue) - The start time is taken from the end\n"
204 "of the previous file. Use it when your continuation run\n"
205 "restarts with t=0 and there is no overlap.\n\n"
206 "l (last) - The time in this file will be changed the\n"
207 "same amount as in the previous. Use it when the time in the\n"
208 "new run continues from the end of the previous one,\n"
209 "since this takes possible overlap into account.\n\n");
211 fprintf(stderr
," File Current start New start\n"
212 "---------------------------------------------------------\n");
214 for(i
=0;i
<nfiles
;i
++) {
215 fprintf(stderr
,"%25s %10.3f ",fnms
[i
],readtime
[i
]);
218 if(NULL
==fgets(inputstring
,STRLEN
-1,stdin
))
220 gmx_fatal(FARGS
,"Error reading user input");
222 inputstring
[strlen(inputstring
)-1]=0;
224 if(inputstring
[0]=='c' || inputstring
[0]=='C') {
225 cont_type
[i
]=TIME_CONTINUE
;
230 else if(inputstring
[0]=='l' ||
231 inputstring
[0]=='L') {
232 cont_type
[i
]=TIME_LAST
;
238 settime
[i
]=strtod(inputstring
,&chptr
);
239 if(chptr
==inputstring
) {
240 fprintf(stderr
,"Try that again: ");
243 cont_type
[i
]=TIME_EXPLICIT
;
249 if(cont_type
[0]!=TIME_EXPLICIT
) {
250 cont_type
[0]=TIME_EXPLICIT
;
255 for(i
=0;i
<nfiles
;i
++)
256 settime
[i
]=readtime
[i
];
258 if(bSort
&& (nfiles
>1))
259 sort_files(fnms
,settime
,nfiles
);
261 fprintf(stderr
,"Sorting disabled.\n");
264 /* Write out the new order and start times */
265 fprintf(stderr
,"\nSummary of files and start times used:\n\n"
267 "-----------------------------------------\n");
268 for(i
=0;i
<nfiles
;i
++)
269 switch(cont_type
[i
]) {
271 fprintf(stderr
,"%25s %10.3f\n",fnms
[i
],settime
[i
]);
274 fprintf(stderr
,"%25s Continue from end of last file\n",fnms
[i
]);
277 fprintf(stderr
,"%25s Change by same amount as last file\n",fnms
[i
]);
280 fprintf(stderr
,"\n");
282 settime
[nfiles
]=FLT_MAX
;
283 cont_type
[nfiles
]=TIME_EXPLICIT
;
284 readtime
[nfiles
]=FLT_MAX
;
288 static void copy_ee(t_energy
*src
, t_energy
*dst
, int nre
)
294 dst
[i
].esum
=src
[i
].esum
;
295 dst
[i
].eav
=src
[i
].eav
;
299 static void update_ee(t_energy
*lastee
,gmx_large_int_t laststep
,
300 t_energy
*startee
,gmx_large_int_t startstep
,
301 t_energy
*ee
, int step
,
302 t_energy
*outee
, int nre
)
305 double sigmacorr
,nom
,denom
;
306 double prestart_esum
;
307 double prestart_sigma
;
311 /* add statistics from earlier file if present */
313 outee
[i
].esum
=lastee
[i
].esum
+ee
[i
].esum
;
314 nom
=(lastee
[i
].esum
*(step
+1)-ee
[i
].esum
*(laststep
));
315 denom
=((step
+1.0)*(laststep
)*(step
+1.0+laststep
));
316 sigmacorr
=nom
*nom
/denom
;
317 outee
[i
].eav
=lastee
[i
].eav
+ee
[i
].eav
+sigmacorr
;
320 /* otherwise just copy to output */
321 outee
[i
].esum
=ee
[i
].esum
;
322 outee
[i
].eav
=ee
[i
].eav
;
325 /* if we didnt start to write at the first frame
326 * we must compensate the statistics for this
327 * there are laststep frames in the earlier file
328 * and step+1 frames in this one.
331 gmx_large_int_t q
=laststep
+step
;
332 gmx_large_int_t p
=startstep
+1;
333 prestart_esum
=startee
[i
].esum
-startee
[i
].e
;
334 sigmacorr
=prestart_esum
-(p
-1)*startee
[i
].e
;
335 prestart_sigma
=startee
[i
].eav
-
336 sigmacorr
*sigmacorr
/(p
*(p
-1));
337 sigmacorr
=prestart_esum
/(p
-1)-
339 outee
[i
].esum
-=prestart_esum
;
341 outee
[i
].eav
=outee
[i
].eav
-prestart_sigma
-
342 sigmacorr
*sigmacorr
*((p
-1)*q
)/(q
-p
+1);
345 if((outee
[i
].eav
/(laststep
+step
+1))<(GMX_REAL_EPS
))
350 static void update_ee_sum(int nre
,
351 gmx_large_int_t
*ee_sum_step
,
352 gmx_large_int_t
*ee_sum_nsteps
,
353 gmx_large_int_t
*ee_sum_nsum
,
355 t_enxframe
*fr
,int out_step
)
357 gmx_large_int_t nsteps
,nsum
,fr_nsum
;
360 nsteps
= *ee_sum_nsteps
;
371 ee_sum
[i
].esum
= fr
->ener
[i
].e
;
376 ee_sum
[i
].esum
= fr
->ener
[i
].esum
;
377 ee_sum
[i
].eav
= fr
->ener
[i
].eav
;
382 } else if (out_step
+ *ee_sum_nsum
- *ee_sum_step
== nsteps
+ fr
->nsteps
) {
386 dsqr(ee_sum
[i
].esum
/nsum
387 - (ee_sum
[i
].esum
+ fr
->ener
[i
].e
)/(nsum
+ 1))*nsum
*(nsum
+ 1);
388 ee_sum
[i
].esum
+= fr
->ener
[i
].e
;
391 for(i
=0; i
<fr
->nre
; i
++) {
394 dsqr(ee_sum
[i
].esum
/nsum
395 - (ee_sum
[i
].esum
+ fr
->ener
[i
].esum
)/(nsum
+ fr
->nsum
))*
396 nsum
*(nsum
+ fr
->nsum
)/(double)fr
->nsum
;
397 ee_sum
[i
].esum
+= fr
->ener
[i
].esum
;
400 nsteps
+= fr
->nsteps
;
404 fprintf(stderr
,"\nWARNING: missing energy sums at time %f\n",fr
->t
);
410 *ee_sum_step
= out_step
;
411 *ee_sum_nsteps
= nsteps
;
415 int gmx_eneconv(int argc
,char *argv
[])
417 const char *desc
[] = {
418 "With [IT]multiple files[it] specified for the [TT]-f[tt] option:[BR]",
419 "Concatenates several energy files in sorted order.",
420 "In the case of double time frames, the one",
421 "in the later file is used. By specifying [TT]-settime[tt] you will be",
422 "asked for the start time of each file. The input files are taken",
423 "from the command line,",
424 "such that the command [TT]eneconv -f *.edr -o fixed.edr[tt] should do",
426 "With [IT]one file[it] specified for [TT]-f[tt]:[BR]",
427 "Reads one energy file and writes another, applying the [TT]-dt[tt],",
428 "[TT]-offset[tt], [TT]-t0[tt] and [TT]-settime[tt] options and",
429 "converting to a different format if necessary (indicated by file",
431 "[TT]-settime[tt] is applied first, then [TT]-dt[tt]/[TT]-offset[tt]",
432 "followed by [TT]-b[tt] and [TT]-e[tt] to select which frames to write."
434 const char *bugs
[] = {
435 "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."
437 ener_file_t in
=NULL
, out
=NULL
;
438 gmx_enxnm_t
*enm
=NULL
;
440 ener_file_t in
,out
=NULL
;
441 gmx_enxnm_t
*enm
=NULL
;
444 gmx_large_int_t ee_sum_step
=0,ee_sum_nsteps
,ee_sum_nsum
;
446 gmx_large_int_t lastfilestep
,laststep
,startstep
,startstep_file
=0;
448 int nre
,nremax
,this_nre
,nfile
,f
,i
,j
,kkk
,nset
,*set
=NULL
;
451 real
*readtime
,*settime
,timestep
,t1
,tadjust
;
452 char inputstring
[STRLEN
],*chptr
,buf
[22],buf2
[22],buf3
[22];
455 gmx_bool bNewFile
,bFirst
,bNewOutput
;
457 gmx_bool warned_about_dh
=FALSE
;
458 t_enxblock
*blocks
=NULL
;
463 { efEDR
, "-f", NULL
, ffRDMULT
},
464 { efEDR
, "-o", "fixed", ffWRITE
},
467 #define NFILE asize(fnm)
469 static real delta_t
=0.0, toffset
=0,scalefac
=1;
470 static gmx_bool bSetTime
=FALSE
;
471 static gmx_bool bSort
=TRUE
,bError
=TRUE
;
472 static real begin
=-1;
474 gmx_bool remove_dh
=FALSE
;
477 { "-b", FALSE
, etREAL
, {&begin
},
478 "First time to use"},
479 { "-e", FALSE
, etREAL
, {&end
},
481 { "-dt", FALSE
, etREAL
, {&delta_t
},
482 "Only write out frame when t MOD dt = offset" },
483 { "-offset", FALSE
, etREAL
, {&toffset
},
484 "Time offset for [TT]-dt[tt] option" },
485 { "-settime", FALSE
, etBOOL
, {&bSetTime
},
486 "Change starting time interactively" },
487 { "-sort", FALSE
, etBOOL
, {&bSort
},
488 "Sort energy files (not frames)"},
489 { "-rmdh", FALSE
, etBOOL
, {&remove_dh
},
490 "Remove free energy block data" },
491 { "-scalefac", FALSE
, etREAL
, {&scalefac
},
492 "Multiply energy component by this factor" },
493 { "-error", FALSE
, etBOOL
, {&bError
},
494 "Stop on errors in the file" }
497 CopyRight(stderr
,argv
[0]);
498 parse_common_args(&argc
,argv
,PCA_BE_NICE
,NFILE
,fnm
,asize(pa
),
499 pa
,asize(desc
),desc
,asize(bugs
),bugs
,&oenv
);
507 laststep
=startstep
=0;
509 nfile
= opt2fns(&fnms
,"-f",NFILE
,fnm
);
512 gmx_fatal(FARGS
,"No input files!");
514 snew(settime
,nfile
+1);
515 snew(readtime
,nfile
+1);
516 snew(cont_type
,nfile
+1);
518 nre
=scan_ene_files(fnms
,nfile
,readtime
,×tep
,&nremax
);
519 edit_files(fnms
,nfile
,readtime
,settime
,cont_type
,bSetTime
,bSort
);
529 snew(fro
->ener
,nremax
);
535 for(f
=0;f
<nfile
;f
++) {
538 in
=open_enx(fnms
[f
],"r");
540 do_enxnms(in
,&this_nre
,&enm
);
543 set
= select_it(nre
,enm
,&nset
);
545 /* write names to the output file */
546 out
=open_enx(opt2fn("-o",NFILE
,fnm
),"w");
547 do_enxnms(out
,&nre
,&enm
);
550 /* start reading from the next file */
551 while ((fro
->t
<= (settime
[f
+1] + GMX_REAL_EPS
)) &&
554 startstep_file
= fr
->step
;
555 tadjust
= settime
[f
] - fr
->t
;
556 if(cont_type
[f
+1]==TIME_LAST
) {
557 settime
[f
+1] = readtime
[f
+1]-readtime
[f
]+settime
[f
];
558 cont_type
[f
+1] = TIME_EXPLICIT
;
563 if (tadjust
+ fr
->t
<= last_t
) {
564 /* Skip this frame, since we already have it / past it */
566 fprintf(debug
,"fr->step %s, fr->t %.4f\n",
567 gmx_step_str(fr
->step
,buf
),fr
->t
);
568 fprintf(debug
,"tadjust %12.6e + fr->t %12.6e <= t %12.6e\n",
569 tadjust
,fr
->t
,last_t
);
574 fro
->step
= lastfilestep
+ fr
->step
- startstep_file
;
575 fro
->t
= tadjust
+ fr
->t
;
577 bWrite
= ((begin
<0 || (begin
>=0 && (fro
->t
>= begin
-GMX_REAL_EPS
))) &&
578 (end
<0 || (end
>=0 && (fro
->t
<= end
+GMX_REAL_EPS
))) &&
579 (fro
->t
<= settime
[f
+1]+0.5*timestep
));
583 "fr->step %s, fr->t %.4f, fro->step %s fro->t %.4f, w %d\n",
584 gmx_step_str(fr
->step
,buf
),fr
->t
,
585 gmx_step_str(fro
->step
,buf2
),fro
->t
,bWrite
);
589 if ((end
> 0) && (fro
->t
> end
+GMX_REAL_EPS
)) {
594 if (fro
->t
>= begin
-GMX_REAL_EPS
) {
597 startstep
= fr
->step
;
600 update_ee_sum(nre
,&ee_sum_step
,&ee_sum_nsteps
,&ee_sum_nsum
,ee_sum
,
605 /* determine if we should write it */
606 if (bWrite
&& (delta_t
==0 || bRmod(fro
->t
,toffset
,delta_t
))) {
607 laststep
= fro
->step
;
611 fprintf(stderr
,"\nContinue writing frames from t=%g, step=%s\n",
612 fro
->t
,gmx_step_str(fro
->step
,buf
));
615 /* Copy the energies */
616 for(i
=0; i
<nre
; i
++) {
617 fro
->ener
[i
].e
= fr
->ener
[i
].e
;
620 fro
->nsteps
= ee_sum_nsteps
;
623 if (ee_sum_nsum
<= 1) {
626 fro
->nsum
= gmx_large_int_to_int(ee_sum_nsum
,
627 "energy average summation");
628 /* Copy the energy sums */
629 for(i
=0; i
<nre
; i
++) {
630 fro
->ener
[i
].esum
= ee_sum
[i
].esum
;
631 fro
->ener
[i
].eav
= ee_sum
[i
].eav
;
634 /* We wrote the energies, so reset the counts */
639 for(kkk
=0; kkk
<nset
; kkk
++) {
640 fro
->ener
[set
[kkk
]].e
*= scalefac
;
642 fro
->ener
[set
[kkk
]].eav
*= scalefac
*scalefac
;
643 fro
->ener
[set
[kkk
]].esum
*= scalefac
;
647 /* Copy restraint stuff */
648 /*fro->ndisre = fr->ndisre;
649 fro->disre_rm3tav = fr->disre_rm3tav;
650 fro->disre_rt = fr->disre_rt;*/
651 fro
->nblock
= fr
->nblock
;
652 /*fro->nr = fr->nr;*/
653 fro
->block
= fr
->block
;
655 /* check if we have blocks with delta_h data and are throwing
662 if (!blocks
|| nblocks_alloc
< fr
->nblock
)
664 /* we pre-allocate the blocks */
665 nblocks_alloc
=fr
->nblock
;
666 snew(blocks
, nblocks_alloc
);
668 nblocks
=0; /* number of blocks so far */
670 for(i
=0;i
<fr
->nblock
;i
++)
672 if ( (fr
->block
[i
].id
!= enxDHCOLL
) &&
673 (fr
->block
[i
].id
!= enxDH
) &&
674 (fr
->block
[i
].id
!= enxDHHIST
) )
676 /* copy everything verbatim */
677 blocks
[nblocks
] = fr
->block
[i
];
681 /* now set the block pointer to the new blocks */
682 fro
->nblock
= nblocks
;
685 else if (delta_t
> 0)
687 if (!warned_about_dh
)
689 for(i
=0;i
<fr
->nblock
;i
++)
691 if (fr
->block
[i
].id
== enxDH
||
692 fr
->block
[i
].id
== enxDHHIST
)
695 if (fr
->block
[i
].id
== enxDH
)
697 size
=fr
->block
[i
].sub
[2].nr
;
705 printf("\nWARNING: %s contains delta H blocks or histograms for which\n"
706 " some data is thrown away on a block-by-block basis, where each block\n"
707 " contains up to %d samples.\n"
708 " This is almost certainly not what you want.\n"
709 " Use the -rmdh option to throw all delta H samples away.\n"
710 " Use g_energy -odh option to extract these samples.\n",
712 warned_about_dh
=TRUE
;
722 if (noutfr
% 1000 == 0)
723 fprintf(stderr
,"Writing frame time %g ",fro
->t
);
729 printf("\nLast step written from %s: t %g, step %s\n",
730 fnms
[f
],last_t
,gmx_step_str(laststep
,buf
));
731 lastfilestep
= laststep
;
733 /* set the next time from the last in previous file */
734 if (cont_type
[f
+1]==TIME_CONTINUE
) {
735 settime
[f
+1] = fro
->t
;
736 /* in this case we have already written the last frame of
737 * previous file, so update begin to avoid doubling it
738 * with the start of the next file
740 begin
= fro
->t
+0.5*timestep
;
741 /* cont_type[f+1]==TIME_EXPLICIT; */
744 if ((fro
->t
< end
) && (f
< nfile
-1) &&
745 (fro
->t
< settime
[f
+1]-1.5*timestep
))
747 "\nWARNING: There might be a gap around t=%g\n",fro
->t
);
749 /* move energies to lastee */
751 free_enxnms(this_nre
,enm
);
753 fprintf(stderr
,"\n");
756 fprintf(stderr
,"No frames written.\n");
758 fprintf(stderr
,"Last frame written was at step %s, time %f\n",
759 gmx_step_str(fro
->step
,buf
),fro
->t
);
760 fprintf(stderr
,"Wrote %d frames\n",noutfr
);