1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * For more info, check our website at http://www.gromacs.org
34 * GROningen Mixture of Alchemy and Childrens' Stories
42 #include "gmx_fatal.h"
48 /* The source code in this file should be thread-safe.
49 Please keep it that way. */
51 /* This number should be increased whenever the file format changes! */
52 static const int enx_version
= 3;
55 /* Stuff for reading pre 4.1 energy files */
57 bool bOldFileOpen
; /* Is this an open old file? */
58 bool bReadFirstStep
; /* Did we read the first step? */
59 int first_step
; /* First step in the energy file */
60 int step_prev
; /* Previous step */
61 int nsum_prev
; /* Previous step sum length */
62 t_energy
*ener_prev
; /* Previous energy sums */
74 void free_enxframe(t_enxframe
*fr
)
81 sfree(fr
->disre_rm3tav
);
84 for(b
=0; b
<fr
->nblock
; b
++)
91 static void edr_strings(XDR
*xdr
,bool bRead
,int file_version
,
92 int n
,gmx_enxnm_t
**nms
)
117 if(!xdr_string(xdr
,&(nm
->name
),STRLEN
))
119 gmx_file("Cannot write energy names to file; maybe you are out of quota?");
121 if (file_version
>= 2)
123 if(!xdr_string(xdr
,&(nm
->unit
),STRLEN
))
125 gmx_file("Cannot write energy names to file; maybe you are out of quota?");
130 nm
->unit
= strdup("kJ/mol");
135 static void gen_units(int n
,char ***units
)
142 (*units
)[i
] = strdup("kJ/mol");
146 void do_enxnms(ener_file_t ef
,int *nre
,gmx_enxnm_t
**nms
)
150 bool bRead
= gmx_fio_getread(ef
->fp
);
154 gmx_fio_select(ef
->fp
);
156 xdr
= gmx_fio_getxdr(ef
->fp
);
158 if (!xdr_int(xdr
,&magic
))
162 gmx_file("Cannot write energy names to file; maybe you are out of quota?");
169 /* Assume this is an old edr format */
172 ef
->eo
.bOldFileOpen
= TRUE
;
173 ef
->eo
.bReadFirstStep
= FALSE
;
174 srenew(ef
->eo
.ener_prev
,*nre
);
178 ef
->eo
.bOldFileOpen
=FALSE
;
182 gmx_fatal(FARGS
,"Energy names magic number mismatch, this is not a GROMACS edr file");
184 file_version
= enx_version
;
185 xdr_int(xdr
,&file_version
);
186 if (file_version
> enx_version
)
188 gmx_fatal(FARGS
,"reading tpx file (%s) version %d with version %d program",gmx_fio_getname(ef
->fp
),file_version
,enx_version
);
192 if (file_version
!= enx_version
)
194 fprintf(stderr
,"Note: enx file_version %d, software version %d\n",
195 file_version
,enx_version
);
198 edr_strings(xdr
,bRead
,file_version
,*nre
,nms
);
201 static bool do_eheader(ener_file_t ef
,int *file_version
,t_enxframe
*fr
,
202 bool bTest
, bool *bOK
)
206 int block
,i
,zero
=0,dum
=0;
207 bool bRead
= gmx_fio_getread(ef
->fp
);
211 /* The original energy frame started with a real,
212 * so we have to use a real for compatibility.
215 if (!do_real(r
)) return FALSE
;
218 /* Assume we are reading an old format */
221 if (!do_int(dum
)) *bOK
= FALSE
;
226 if (!do_int (magic
)) *bOK
= FALSE
;
227 if (magic
!= -7777777)
229 gmx_fatal(FARGS
,"Energy header magic number mismatch, this is not a GROMACS edr file");
231 *file_version
= enx_version
;
232 if (!do_int (*file_version
)) *bOK
= FALSE
;
233 if (*bOK
&& *file_version
> enx_version
)
235 gmx_fatal(FARGS
,"reading tpx file (%s) version %d with version %d program",gmx_fio_getname(ef
->fp
),file_version
,enx_version
);
237 if (!do_double(fr
->t
)) *bOK
= FALSE
;
238 if (!do_gmx_large_int(fr
->step
)) *bOK
= FALSE
;
239 if (!bRead
&& fr
->nsum
== 1) {
240 /* Do not store sums of length 1,
241 * since this does not add information.
243 if (!do_int (zero
)) *bOK
= FALSE
;
245 if (!do_int (fr
->nsum
)) *bOK
= FALSE
;
247 if (*file_version
>= 3)
249 do_gmx_large_int(fr
->nsteps
);
253 fr
->nsteps
= max(1,fr
->nsum
);
256 if (!do_int (fr
->nre
)) *bOK
= FALSE
;
257 if (!do_int (fr
->ndisre
)) *bOK
= FALSE
;
258 if (!do_int (fr
->nblock
)) *bOK
= FALSE
;
260 if (*bOK
&& bRead
&& fr
->nblock
>fr
->nr_alloc
)
262 srenew(fr
->nr
,fr
->nblock
);
263 srenew(fr
->b_alloc
,fr
->nblock
);
264 srenew(fr
->block
,fr
->nblock
);
265 for(i
=fr
->nr_alloc
; i
<fr
->nblock
; i
++) {
269 fr
->nr_alloc
= fr
->nblock
;
271 for(block
=0; block
<fr
->nblock
; block
++)
273 if (!do_int (fr
->nr
[block
]))
278 if (!do_int (fr
->e_size
)) *bOK
= FALSE
;
279 if (!do_int (fr
->d_size
)) *bOK
= FALSE
;
280 /* Do a dummy int to keep the format compatible with the old code */
281 if (!do_int (dum
)) *bOK
= FALSE
;
284 if (*bOK
&& *file_version
== 1 && !bTest
)
287 if (fp
>= ener_old_nalloc
)
289 gmx_incons("Problem with reading old format energy files");
293 if (!ef
->eo
.bReadFirstStep
)
295 ef
->eo
.bReadFirstStep
= TRUE
;
296 ef
->eo
.first_step
= fr
->step
;
297 ef
->eo
.nsum_prev
= 0;
300 fr
->nsum
= fr
->step
- ef
->eo
.first_step
+ 1;
301 fr
->nsteps
= fr
->nsum
;
307 void free_enxnms(int n
,gmx_enxnm_t
*nms
)
320 void close_enx(ener_file_t ef
)
322 if(gmx_fio_close(ef
->fp
) != 0)
324 gmx_file("Cannot close energy file; it might be corrupt, or maybe you are out of quota?");
328 static bool empty_file(const char *fn
)
335 fp
= gmx_fio_fopen(fn
,"r");
336 ret
= fread(&dum
,sizeof(dum
),1,fp
);
344 ener_file_t
open_enx(const char *fn
,const char *mode
)
347 gmx_enxnm_t
*nms
=NULL
;
351 struct ener_file
*ef
;
356 ef
->fp
=gmx_fio_open(fn
,mode
);
357 gmx_fio_select(ef
->fp
);
358 gmx_fio_setprecision(ef
->fp
,FALSE
);
359 do_enxnms(ef
,&nre
,&nms
);
361 do_eheader(ef
,&file_version
,fr
,TRUE
,&bDum
);
364 gmx_file("Cannot read energy file header. Corrupt file?");
367 /* Now check whether this file is in single precision */
368 if (((fr
->e_size
&& (fr
->nre
== nre
) &&
369 (nre
*4*(long int)sizeof(float) == fr
->e_size
)) ||
371 (fr
->ndisre
*(long int)sizeof(float)*2+(long int)sizeof(int)
373 fprintf(stderr
,"Opened %s as single precision energy file\n",fn
);
374 free_enxnms(nre
,nms
);
377 gmx_fio_rewind(ef
->fp
);
378 gmx_fio_select(ef
->fp
);
379 gmx_fio_setprecision(ef
->fp
,TRUE
);
380 do_enxnms(ef
,&nre
,&nms
);
381 do_eheader(ef
,&file_version
,fr
,TRUE
,&bDum
);
384 gmx_file("Cannot write energy file header; maybe you are out of quota?");
387 if (((fr
->e_size
&& (fr
->nre
== nre
) &&
388 (nre
*4*(long int)sizeof(double) == fr
->e_size
)) ||
390 (fr
->ndisre
*(long int)sizeof(double)*2+
391 (long int)sizeof(int) ==
393 fprintf(stderr
,"Opened %s as double precision energy file\n",
397 gmx_fatal(FARGS
,"File %s is empty",fn
);
399 gmx_fatal(FARGS
,"Energy file %s not recognized, maybe different CPU?",
402 free_enxnms(nre
,nms
);
406 gmx_fio_rewind(ef
->fp
);
409 ef
->fp
= gmx_fio_open(fn
,mode
);
417 int enx_file_pointer(const ener_file_t ef
)
422 static void convert_full_sums(ener_old_t
*ener_old
,t_enxframe
*fr
)
426 double esum_all
,eav_all
;
432 for(i
=0; i
<fr
->nre
; i
++)
434 if (fr
->ener
[i
].e
!= 0) ne
++;
435 if (fr
->ener
[i
].esum
!= 0) ns
++;
437 if (ne
> 0 && ns
== 0)
439 /* We do not have all energy sums */
444 /* Convert old full simulation sums to sums between energy frames */
445 nstep_all
= fr
->step
- ener_old
->first_step
+ 1;
446 if (fr
->nsum
> 1 && fr
->nsum
== nstep_all
&& ener_old
->nsum_prev
> 0)
448 /* Set the new sum length: the frame step difference */
449 fr
->nsum
= fr
->step
- ener_old
->step_prev
;
450 for(i
=0; i
<fr
->nre
; i
++)
452 esum_all
= fr
->ener
[i
].esum
;
453 eav_all
= fr
->ener
[i
].eav
;
454 fr
->ener
[i
].esum
= esum_all
- ener_old
->ener_prev
[i
].esum
;
455 fr
->ener
[i
].eav
= eav_all
- ener_old
->ener_prev
[i
].eav
456 - dsqr(ener_old
->ener_prev
[i
].esum
/(nstep_all
- fr
->nsum
)
457 - esum_all
/nstep_all
)*
458 (nstep_all
- fr
->nsum
)*nstep_all
/(double)fr
->nsum
;
459 ener_old
->ener_prev
[i
].esum
= esum_all
;
460 ener_old
->ener_prev
[i
].eav
= eav_all
;
462 ener_old
->nsum_prev
= nstep_all
;
464 else if (fr
->nsum
> 0)
466 if (fr
->nsum
!= nstep_all
)
468 fprintf(stderr
,"\nWARNING: something is wrong with the energy sums, will not use exact averages\n");
469 ener_old
->nsum_prev
= 0;
473 ener_old
->nsum_prev
= nstep_all
;
475 /* Copy all sums to ener_prev */
476 for(i
=0; i
<fr
->nre
; i
++)
478 ener_old
->ener_prev
[i
].esum
= fr
->ener
[i
].esum
;
479 ener_old
->ener_prev
[i
].eav
= fr
->ener
[i
].eav
;
483 ener_old
->step_prev
= fr
->step
;
486 bool do_enx(ener_file_t ef
,t_enxframe
*fr
)
490 bool bRead
,bOK
,bOK1
,bSane
;
495 bRead
= gmx_fio_getread(ef
->fp
);
498 fr
->e_size
= fr
->nre
*sizeof(fr
->ener
[0].e
)*4;
499 fr
->d_size
= fr
->ndisre
*(sizeof(fr
->disre_rm3tav
[0]) +
500 sizeof(fr
->disre_rt
[0]));
502 gmx_fio_select(ef
->fp
);
504 if (!do_eheader(ef
,&file_version
,fr
,FALSE
,&bOK
))
508 fprintf(stderr
,"\rLast energy frame read %d time %8.3f ",
509 ef
->framenr
-1,ef
->frametime
);
513 "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
519 gmx_file("Cannot write energy file header; maybe you are out of quota?");
525 if ((ef
->framenr
< 20 || ef
->framenr
% 10 == 0) &&
526 (ef
->framenr
< 200 || ef
->framenr
% 100 == 0) &&
527 (ef
->framenr
< 2000 || ef
->framenr
% 1000 == 0))
529 fprintf(stderr
,"\rReading energy frame %6d time %8.3f ",
533 ef
->frametime
= fr
->t
;
535 /* Check sanity of this header */
536 bSane
= (fr
->nre
> 0 || fr
->ndisre
> 0);
537 for(block
=0; block
<fr
->nblock
; block
++)
539 bSane
= bSane
|| (fr
->nr
[block
] > 0);
541 if (!((fr
->step
>= 0) && bSane
))
543 fprintf(stderr
,"\nWARNING: there may be something wrong with energy file %s\n",
544 gmx_fio_getname(ef
->fp
));
545 fprintf(stderr
,"Found: step=%s, nre=%d, ndisre=%d, nblock=%d, time=%g.\n"
546 "Trying to skip frame expect a crash though\n",
547 gmx_step_str(fr
->step
,buf
),fr
->nre
,fr
->ndisre
,fr
->nblock
,fr
->t
);
549 if (bRead
&& fr
->nre
> fr
->e_alloc
)
551 srenew(fr
->ener
,fr
->nre
);
552 for(i
=fr
->e_alloc
; (i
<fr
->nre
); i
++)
556 fr
->ener
[i
].esum
= 0;
558 fr
->e_alloc
= fr
->nre
;
561 for(i
=0; i
<fr
->nre
; i
++)
563 bOK
= bOK
&& do_real(fr
->ener
[i
].e
);
565 /* Do not store sums of length 1,
566 * since this does not add information.
568 if (file_version
== 1 ||
569 (bRead
&& fr
->nsum
> 0) || fr
->nsum
> 1)
571 tmp1
= fr
->ener
[i
].eav
;
572 bOK
= bOK
&& do_real(tmp1
);
574 fr
->ener
[i
].eav
= tmp1
;
576 /* This is to save only in single precision (unless compiled in DP) */
577 tmp2
= fr
->ener
[i
].esum
;
578 bOK
= bOK
&& do_real(tmp2
);
580 fr
->ener
[i
].esum
= tmp2
;
582 if (file_version
== 1)
584 /* Old, unused real */
586 bOK
= bOK
&& do_real(rdum
);
591 /* Here we can not check for file_version==1, since one could have
592 * continued an old format simulation with a new one with mdrun -append.
594 if (bRead
&& ef
->eo
.bOldFileOpen
)
596 /* Convert old full simulation sums to sums between energy frames */
597 convert_full_sums(&(ef
->eo
),fr
);
601 if (bRead
&& fr
->ndisre
>fr
->d_alloc
)
603 srenew(fr
->disre_rm3tav
,fr
->ndisre
);
604 srenew(fr
->disre_rt
,fr
->ndisre
);
605 fr
->d_alloc
= fr
->ndisre
;
607 ndo_real(fr
->disre_rm3tav
,fr
->ndisre
,bOK1
);
609 ndo_real(fr
->disre_rt
,fr
->ndisre
,bOK1
);
612 for(block
=0; block
<fr
->nblock
; block
++)
614 if (bRead
&& fr
->nr
[block
]>fr
->b_alloc
[block
])
616 srenew(fr
->block
[block
],fr
->nr
[block
]);
617 fr
->b_alloc
[block
] = fr
->nr
[block
];
619 ndo_real(fr
->block
[block
],fr
->nr
[block
],bOK1
);
625 if( gmx_fio_flush(ef
->fp
) != 0)
627 gmx_file("Cannot write energy file; maybe you are out of quota?");
635 fprintf(stderr
,"\nLast energy frame read %d",
637 fprintf(stderr
,"\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
642 gmx_fatal(FARGS
,"could not write energies");
650 static real
find_energy(const char *name
, int nre
, gmx_enxnm_t
*enm
,
657 if (strcmp(enm
[i
].name
,name
) == 0)
659 return fr
->ener
[i
].e
;
663 gmx_fatal(FARGS
,"Could not find energy term named '%s'",name
);
669 void get_enx_state(const char *fn
, real t
, gmx_groups_t
*groups
, t_inputrec
*ir
,
672 /* Should match the names in mdebin.c */
673 static const char *boxvel_nm
[] = {
674 "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
675 "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
678 static const char *pcouplmu_nm
[] = {
679 "Pcoupl-Mu-XX", "Pcoupl-Mu-YY", "Pcoupl-Mu-ZZ",
680 "Pcoupl-Mu-YX", "Pcoupl-Mu-ZX", "Pcoupl-Mu-ZY"
682 int ind0
[] = { XX
,YY
,ZZ
,YY
,ZZ
,ZZ
};
683 int ind1
[] = { XX
,YY
,ZZ
,XX
,XX
,YY
};
685 int nre
,nfr
,i
,ni
,npcoupl
;
691 in
= open_enx(fn
,"r");
692 do_enxnms(in
,&nre
,&enm
);
695 while ((nfr
==0 || fr
->t
!= t
) && do_enx(in
,fr
)) {
699 fprintf(stderr
,"\n");
701 if (nfr
== 0 || fr
->t
!= t
)
702 gmx_fatal(FARGS
,"Could not find frame with time %f in '%s'",t
,fn
);
704 npcoupl
= TRICLINIC(ir
->compress
) ? 6 : 3;
705 if (ir
->epc
== epcPARRINELLORAHMAN
) {
706 clear_mat(state
->boxv
);
707 for(i
=0; i
<npcoupl
; i
++) {
708 state
->boxv
[ind0
[i
]][ind1
[i
]] =
709 find_energy(boxvel_nm
[i
],nre
,enm
,fr
);
711 fprintf(stderr
,"\nREAD %d BOX VELOCITIES FROM %s\n\n",npcoupl
,fn
);
714 if (ir
->etc
== etcNOSEHOOVER
) {
715 for(i
=0; i
<state
->ngtc
; i
++) {
716 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
717 sprintf(buf
,"Xi-%s",*(groups
->grpname
[ni
]));
718 state
->nosehoover_xi
[i
] = find_energy(buf
,nre
,enm
,fr
);
720 fprintf(stderr
,"\nREAD %d NOSE-HOOVER Xi's FROM %s\n\n",state
->ngtc
,fn
);
723 free_enxnms(nre
,enm
);