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"
49 /* The source code in this file should be thread-safe.
50 Please keep it that way. */
52 /* This number should be increased whenever the file format changes! */
53 static const int enx_version
= 5;
55 const char *enx_block_id_name
[] = {
56 "Averaged orientation restraints",
57 "Instantaneous orientation restraints",
58 "Orientation restraint order tensor(s)",
59 "Distance restraints",
66 /* Stuff for reading pre 4.1 energy files */
68 gmx_bool bOldFileOpen
; /* Is this an open old file? */
69 gmx_bool bReadFirstStep
; /* Did we read the first step? */
70 int first_step
; /* First step in the energy file */
71 int step_prev
; /* Previous step */
72 int nsum_prev
; /* Previous step sum length */
73 t_energy
*ener_prev
; /* Previous energy sums */
84 static void enxsubblock_init(t_enxsubblock
*sb
)
88 sb
->type
=xdr_datatype_double
;
90 sb
->type
=xdr_datatype_float
;
106 static void enxsubblock_free(t_enxsubblock
*sb
)
142 for(i
=0;i
<sb
->sval_alloc
;i
++)
155 /* allocate the appropriate amount of memory for the given type and nr */
156 static void enxsubblock_alloc(t_enxsubblock
*sb
)
158 /* allocate the appropriate amount of memory */
161 case xdr_datatype_float
:
162 if (sb
->nr
> sb
->fval_alloc
)
164 srenew(sb
->fval
, sb
->nr
);
165 sb
->fval_alloc
=sb
->nr
;
168 case xdr_datatype_double
:
169 if (sb
->nr
> sb
->dval_alloc
)
171 srenew(sb
->dval
, sb
->nr
);
172 sb
->dval_alloc
=sb
->nr
;
175 case xdr_datatype_int
:
176 if (sb
->nr
> sb
->ival_alloc
)
178 srenew(sb
->ival
, sb
->nr
);
179 sb
->ival_alloc
=sb
->nr
;
182 case xdr_datatype_large_int
:
183 if (sb
->nr
> sb
->lval_alloc
)
185 srenew(sb
->lval
, sb
->nr
);
186 sb
->lval_alloc
=sb
->nr
;
189 case xdr_datatype_char
:
190 if (sb
->nr
> sb
->cval_alloc
)
192 srenew(sb
->cval
, sb
->nr
);
193 sb
->cval_alloc
=sb
->nr
;
196 case xdr_datatype_string
:
197 if (sb
->nr
> sb
->sval_alloc
)
201 srenew(sb
->sval
, sb
->nr
);
202 for(i
=sb
->sval_alloc
;i
<sb
->nr
;i
++)
206 sb
->sval_alloc
=sb
->nr
;
210 gmx_incons("Unknown block type: this file is corrupted or from the future");
214 static void enxblock_init(t_enxblock
*eb
)
222 static void enxblock_free(t_enxblock
*eb
)
224 if (eb
->nsub_alloc
>0)
227 for(i
=0;i
<eb
->nsub_alloc
;i
++)
229 enxsubblock_free(&(eb
->sub
[i
]));
237 void init_enxframe(t_enxframe
*fr
)
253 void free_enxframe(t_enxframe
*fr
)
261 for(b
=0; b
<fr
->nblock_alloc
; b
++)
263 enxblock_free(&(fr
->block
[b
]));
268 void add_blocks_enxframe(t_enxframe
*fr
, int n
)
271 if (n
> fr
->nblock_alloc
)
275 srenew(fr
->block
, n
);
276 for(b
=fr
->nblock_alloc
;b
<fr
->nblock
;b
++)
278 enxblock_init(&(fr
->block
[b
]));
284 t_enxblock
*find_block_id_enxframe(t_enxframe
*ef
, int id
, t_enxblock
*prev
)
291 starti
=(prev
- ef
->block
) + 1;
293 for(i
=starti
; i
<ef
->nblock
; i
++)
295 if (ef
->block
[i
].id
== id
)
296 return &(ef
->block
[i
]);
301 void add_subblocks_enxblock(t_enxblock
*eb
, int n
)
304 if (eb
->nsub
> eb
->nsub_alloc
)
309 for(b
=eb
->nsub_alloc
; b
<n
; b
++)
311 enxsubblock_init(&(eb
->sub
[b
]));
317 static void enx_warning(const char *msg
)
319 if (getenv("GMX_ENX_NO_FATAL") != NULL
)
325 gmx_fatal(FARGS
,"%s\n%s",
327 "If you want to use the correct frames before the corrupted frame and avoid this fatal error set the env.var. GMX_ENX_NO_FATAL");
331 static void edr_strings(XDR
*xdr
,gmx_bool bRead
,int file_version
,
332 int n
,gmx_enxnm_t
**nms
)
357 if(!xdr_string(xdr
,&(nm
->name
),STRLEN
))
359 gmx_file("Cannot write energy names to file; maybe you are out of quota?");
361 if (file_version
>= 2)
363 if(!xdr_string(xdr
,&(nm
->unit
),STRLEN
))
365 gmx_file("Cannot write energy names to file; maybe you are out of quota?");
370 nm
->unit
= strdup("kJ/mol");
375 static void gen_units(int n
,char ***units
)
382 (*units
)[i
] = strdup("kJ/mol");
386 void do_enxnms(ener_file_t ef
,int *nre
,gmx_enxnm_t
**nms
)
390 gmx_bool bRead
= gmx_fio_getread(ef
->fio
);
394 gmx_fio_checktype(ef
->fio
);
396 xdr
= gmx_fio_getxdr(ef
->fio
);
398 if (!xdr_int(xdr
,&magic
))
402 gmx_file("Cannot write energy names to file; maybe you are out of quota?");
409 /* Assume this is an old edr format */
412 ef
->eo
.bOldFileOpen
= TRUE
;
413 ef
->eo
.bReadFirstStep
= FALSE
;
414 srenew(ef
->eo
.ener_prev
,*nre
);
418 ef
->eo
.bOldFileOpen
=FALSE
;
422 gmx_fatal(FARGS
,"Energy names magic number mismatch, this is not a GROMACS edr file");
424 file_version
= enx_version
;
425 xdr_int(xdr
,&file_version
);
426 if (file_version
> enx_version
)
428 gmx_fatal(FARGS
,"reading tpx file (%s) version %d with version %d program",gmx_fio_getname(ef
->fio
),file_version
,enx_version
);
432 if (file_version
!= enx_version
)
434 fprintf(stderr
,"Note: enx file_version %d, software version %d\n",
435 file_version
,enx_version
);
438 edr_strings(xdr
,bRead
,file_version
,*nre
,nms
);
441 static gmx_bool
do_eheader(ener_file_t ef
,int *file_version
,t_enxframe
*fr
,
442 int nre_test
,gmx_bool
*bWrongPrecision
,gmx_bool
*bOK
)
445 real first_real_to_check
;
446 int b
,i
,zero
=0,dum
=0;
447 gmx_bool bRead
= gmx_fio_getread(ef
->fio
);
452 xdr_datatype dtreal
=xdr_datatype_float
;
454 xdr_datatype dtreal
=xdr_datatype_double
;
459 *bWrongPrecision
= FALSE
;
463 /* The original energy frame started with a real,
464 * so we have to use a real for compatibility.
465 * This is VERY DIRTY code, since do_eheader can be called
466 * with the wrong precision set and then we could read r > -1e10,
467 * while actually the intention was r < -1e10.
468 * When nre_test >= 0, do_eheader should therefore terminate
469 * before the number of i/o calls starts depending on what has been read
470 * (which is the case for for instance the block sizes for variable
471 * number of blocks, where this number is read before).
473 first_real_to_check
= -2e10
;
474 if (!gmx_fio_do_real(ef
->fio
, first_real_to_check
))
478 if (first_real_to_check
> -1e10
)
480 /* Assume we are reading an old format */
482 fr
->t
= first_real_to_check
;
483 if (!gmx_fio_do_int(ef
->fio
, dum
)) *bOK
= FALSE
;
488 if (!gmx_fio_do_int(ef
->fio
, magic
)) *bOK
= FALSE
;
489 if (magic
!= -7777777)
491 enx_warning("Energy header magic number mismatch, this is not a GROMACS edr file");
495 *file_version
= enx_version
;
496 if (!gmx_fio_do_int(ef
->fio
, *file_version
)) *bOK
= FALSE
;
497 if (*bOK
&& *file_version
> enx_version
)
499 gmx_fatal(FARGS
,"reading tpx file (%s) version %d with version %d program",gmx_fio_getname(ef
->fio
),file_version
,enx_version
);
501 if (!gmx_fio_do_double(ef
->fio
, fr
->t
)) *bOK
= FALSE
;
502 if (!gmx_fio_do_gmx_large_int(ef
->fio
, fr
->step
)) *bOK
= FALSE
;
503 if (!bRead
&& fr
->nsum
== 1) {
504 /* Do not store sums of length 1,
505 * since this does not add information.
507 if (!gmx_fio_do_int(ef
->fio
, zero
)) *bOK
= FALSE
;
509 if (!gmx_fio_do_int(ef
->fio
, fr
->nsum
)) *bOK
= FALSE
;
511 if (*file_version
>= 3)
513 if (!gmx_fio_do_gmx_large_int(ef
->fio
, fr
->nsteps
)) *bOK
= FALSE
;
517 fr
->nsteps
= max(1,fr
->nsum
);
519 if (*file_version
>= 5)
521 if (!gmx_fio_do_double(ef
->fio
, fr
->dt
)) *bOK
= FALSE
;
528 if (!gmx_fio_do_int(ef
->fio
, fr
->nre
)) *bOK
= FALSE
;
529 if (*file_version
< 4)
531 if (!gmx_fio_do_int(ef
->fio
, ndisre
)) *bOK
= FALSE
;
535 /* now reserved for possible future use */
536 if (!gmx_fio_do_int(ef
->fio
, dum
)) *bOK
= FALSE
;
539 if (!gmx_fio_do_int(ef
->fio
, fr
->nblock
)) *bOK
= FALSE
;
540 if (fr
->nblock
< 0) *bOK
=FALSE
;
544 if (*file_version
>= 4)
546 enx_warning("Distance restraint blocks in old style in new style file");
554 /* Frames could have nre=0, so we can not rely only on the fr->nre check */
555 if (bRead
&& nre_test
>= 0 &&
556 ((fr
->nre
> 0 && fr
->nre
!= nre_test
) ||
557 fr
->nre
< 0 || ndisre
< 0 || fr
->nblock
< 0))
559 *bWrongPrecision
= TRUE
;
563 /* we now know what these should be, or we've already bailed out because
564 of wrong precision */
565 if ( *file_version
==1 && (fr
->t
< 0 || fr
->t
> 1e20
|| fr
->step
< 0 ) )
567 enx_warning("edr file with negative step number or unreasonable time (and without version number).");
575 add_blocks_enxframe(fr
, fr
->nblock
);
581 /* sub[0] is the instantaneous data, sub[1] is time averaged */
582 add_subblocks_enxblock(&(fr
->block
[0]), 2);
583 fr
->block
[0].id
=enxDISRE
;
584 fr
->block
[0].sub
[0].nr
=ndisre
;
585 fr
->block
[0].sub
[1].nr
=ndisre
;
586 fr
->block
[0].sub
[0].type
=dtreal
;
587 fr
->block
[0].sub
[1].type
=dtreal
;
591 /* read block header info */
592 for(b
=startb
; b
<fr
->nblock
; b
++)
596 /* blocks in old version files always have 1 subblock that
597 consists of reals. */
602 add_subblocks_enxblock(&(fr
->block
[b
]), 1);
606 if (fr
->block
[b
].nsub
!= 1)
608 gmx_incons("Writing an old version .edr file with too many subblocks");
610 if (fr
->block
[b
].sub
[0].type
!= dtreal
)
612 gmx_incons("Writing an old version .edr file the wrong subblock type");
615 nrint
= fr
->block
[b
].sub
[0].nr
;
617 if (!gmx_fio_do_int(ef
->fio
, nrint
))
621 fr
->block
[b
].id
= b
- startb
;
622 fr
->block
[b
].sub
[0].nr
= nrint
;
623 fr
->block
[b
].sub
[0].type
= dtreal
;
628 /* in the new version files, the block header only contains
629 the ID and the number of subblocks */
630 int nsub
=fr
->block
[b
].nsub
;
631 *bOK
= *bOK
&& gmx_fio_do_int(ef
->fio
, fr
->block
[b
].id
);
632 *bOK
= *bOK
&& gmx_fio_do_int(ef
->fio
, nsub
);
634 fr
->block
[b
].nsub
=nsub
;
636 add_subblocks_enxblock(&(fr
->block
[b
]), nsub
);
638 /* read/write type & size for each subblock */
641 t_enxsubblock
*sub
=&(fr
->block
[b
].sub
[i
]); /* shortcut */
642 int typenr
=sub
->type
;
644 *bOK
=*bOK
&& gmx_fio_do_int(ef
->fio
, typenr
);
645 *bOK
=*bOK
&& gmx_fio_do_int(ef
->fio
, sub
->nr
);
647 sub
->type
= (xdr_datatype
)typenr
;
651 if (!gmx_fio_do_int(ef
->fio
, fr
->e_size
)) *bOK
= FALSE
;
653 /* now reserved for possible future use */
654 if (!gmx_fio_do_int(ef
->fio
, dum
)) *bOK
= FALSE
;
656 /* Do a dummy int to keep the format compatible with the old code */
657 if (!gmx_fio_do_int(ef
->fio
, dum
)) *bOK
= FALSE
;
659 if (*bOK
&& *file_version
== 1 && nre_test
< 0)
662 if (fp
>= ener_old_nalloc
)
664 gmx_incons("Problem with reading old format energy files");
668 if (!ef
->eo
.bReadFirstStep
)
670 ef
->eo
.bReadFirstStep
= TRUE
;
671 ef
->eo
.first_step
= fr
->step
;
672 ef
->eo
.step_prev
= fr
->step
;
673 ef
->eo
.nsum_prev
= 0;
676 fr
->nsum
= fr
->step
- ef
->eo
.first_step
+ 1;
677 fr
->nsteps
= fr
->step
- ef
->eo
.step_prev
;
684 void free_enxnms(int n
,gmx_enxnm_t
*nms
)
697 void close_enx(ener_file_t ef
)
699 if(gmx_fio_close(ef
->fio
) != 0)
701 gmx_file("Cannot close energy file; it might be corrupt, or maybe you are out of quota?");
705 static gmx_bool
empty_file(const char *fn
)
712 fp
= gmx_fio_fopen(fn
,"r");
713 ret
= fread(&dum
,sizeof(dum
),1,fp
);
721 ener_file_t
open_enx(const char *fn
,const char *mode
)
724 gmx_enxnm_t
*nms
=NULL
;
727 gmx_bool bWrongPrecision
,bOK
=TRUE
;
728 struct ener_file
*ef
;
733 ef
->fio
=gmx_fio_open(fn
,mode
);
734 gmx_fio_checktype(ef
->fio
);
735 gmx_fio_setprecision(ef
->fio
,FALSE
);
736 do_enxnms(ef
,&nre
,&nms
);
738 do_eheader(ef
,&file_version
,fr
,nre
,&bWrongPrecision
,&bOK
);
741 gmx_file("Cannot read energy file header. Corrupt file?");
744 /* Now check whether this file is in single precision */
745 if (!bWrongPrecision
&&
746 ((fr
->e_size
&& (fr
->nre
== nre
) &&
747 (nre
*4*(long int)sizeof(float) == fr
->e_size
)) ) )
749 fprintf(stderr
,"Opened %s as single precision energy file\n",fn
);
750 free_enxnms(nre
,nms
);
754 gmx_fio_rewind(ef
->fio
);
755 gmx_fio_checktype(ef
->fio
);
756 gmx_fio_setprecision(ef
->fio
,TRUE
);
757 do_enxnms(ef
,&nre
,&nms
);
758 do_eheader(ef
,&file_version
,fr
,nre
,&bWrongPrecision
,&bOK
);
761 gmx_file("Cannot write energy file header; maybe you are out of quota?");
764 if (((fr
->e_size
&& (fr
->nre
== nre
) &&
765 (nre
*4*(long int)sizeof(double) == fr
->e_size
)) ))
766 fprintf(stderr
,"Opened %s as double precision energy file\n",
770 gmx_fatal(FARGS
,"File %s is empty",fn
);
772 gmx_fatal(FARGS
,"Energy file %s not recognized, maybe different CPU?",
775 free_enxnms(nre
,nms
);
779 gmx_fio_rewind(ef
->fio
);
782 ef
->fio
= gmx_fio_open(fn
,mode
);
789 t_fileio
*enx_file_pointer(const ener_file_t ef
)
794 static void convert_full_sums(ener_old_t
*ener_old
,t_enxframe
*fr
)
798 double esum_all
,eav_all
;
804 for(i
=0; i
<fr
->nre
; i
++)
806 if (fr
->ener
[i
].e
!= 0) ne
++;
807 if (fr
->ener
[i
].esum
!= 0) ns
++;
809 if (ne
> 0 && ns
== 0)
811 /* We do not have all energy sums */
816 /* Convert old full simulation sums to sums between energy frames */
817 nstep_all
= fr
->step
- ener_old
->first_step
+ 1;
818 if (fr
->nsum
> 1 && fr
->nsum
== nstep_all
&& ener_old
->nsum_prev
> 0)
820 /* Set the new sum length: the frame step difference */
821 fr
->nsum
= fr
->step
- ener_old
->step_prev
;
822 for(i
=0; i
<fr
->nre
; i
++)
824 esum_all
= fr
->ener
[i
].esum
;
825 eav_all
= fr
->ener
[i
].eav
;
826 fr
->ener
[i
].esum
= esum_all
- ener_old
->ener_prev
[i
].esum
;
827 fr
->ener
[i
].eav
= eav_all
- ener_old
->ener_prev
[i
].eav
828 - dsqr(ener_old
->ener_prev
[i
].esum
/(nstep_all
- fr
->nsum
)
829 - esum_all
/nstep_all
)*
830 (nstep_all
- fr
->nsum
)*nstep_all
/(double)fr
->nsum
;
831 ener_old
->ener_prev
[i
].esum
= esum_all
;
832 ener_old
->ener_prev
[i
].eav
= eav_all
;
834 ener_old
->nsum_prev
= nstep_all
;
836 else if (fr
->nsum
> 0)
838 if (fr
->nsum
!= nstep_all
)
840 fprintf(stderr
,"\nWARNING: something is wrong with the energy sums, will not use exact averages\n");
841 ener_old
->nsum_prev
= 0;
845 ener_old
->nsum_prev
= nstep_all
;
847 /* Copy all sums to ener_prev */
848 for(i
=0; i
<fr
->nre
; i
++)
850 ener_old
->ener_prev
[i
].esum
= fr
->ener
[i
].esum
;
851 ener_old
->ener_prev
[i
].eav
= fr
->ener
[i
].eav
;
855 ener_old
->step_prev
= fr
->step
;
858 gmx_bool
do_enx(ener_file_t ef
,t_enxframe
*fr
)
862 gmx_bool bRead
,bOK
,bOK1
,bSane
;
868 bRead
= gmx_fio_getread(ef
->fio
);
871 fr
->e_size
= fr
->nre
*sizeof(fr
->ener
[0].e
)*4;
872 /*d_size = fr->ndisre*(sizeof(real)*2);*/
874 gmx_fio_checktype(ef
->fio
);
876 if (!do_eheader(ef
,&file_version
,fr
,-1,NULL
,&bOK
))
880 fprintf(stderr
,"\rLast energy frame read %d time %8.3f ",
881 ef
->framenr
-1,ef
->frametime
);
885 "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
891 gmx_file("Cannot write energy file header; maybe you are out of quota?");
897 if ((ef
->framenr
< 20 || ef
->framenr
% 10 == 0) &&
898 (ef
->framenr
< 200 || ef
->framenr
% 100 == 0) &&
899 (ef
->framenr
< 2000 || ef
->framenr
% 1000 == 0))
901 fprintf(stderr
,"\rReading energy frame %6d time %8.3f ",
905 ef
->frametime
= fr
->t
;
907 /* Check sanity of this header */
908 bSane
= fr
->nre
> 0 ;
909 for(b
=0; b
<fr
->nblock
; b
++)
911 bSane
= bSane
|| (fr
->block
[b
].nsub
> 0);
913 if (!((fr
->step
>= 0) && bSane
))
915 fprintf(stderr
,"\nWARNING: there may be something wrong with energy file %s\n",
916 gmx_fio_getname(ef
->fio
));
917 fprintf(stderr
,"Found: step=%s, nre=%d, nblock=%d, time=%g.\n"
918 "Trying to skip frame expect a crash though\n",
919 gmx_step_str(fr
->step
,buf
),fr
->nre
,fr
->nblock
,fr
->t
);
921 if (bRead
&& fr
->nre
> fr
->e_alloc
)
923 srenew(fr
->ener
,fr
->nre
);
924 for(i
=fr
->e_alloc
; (i
<fr
->nre
); i
++)
928 fr
->ener
[i
].esum
= 0;
930 fr
->e_alloc
= fr
->nre
;
933 for(i
=0; i
<fr
->nre
; i
++)
935 bOK
= bOK
&& gmx_fio_do_real(ef
->fio
, fr
->ener
[i
].e
);
937 /* Do not store sums of length 1,
938 * since this does not add information.
940 if (file_version
== 1 ||
941 (bRead
&& fr
->nsum
> 0) || fr
->nsum
> 1)
943 tmp1
= fr
->ener
[i
].eav
;
944 bOK
= bOK
&& gmx_fio_do_real(ef
->fio
, tmp1
);
946 fr
->ener
[i
].eav
= tmp1
;
948 /* This is to save only in single precision (unless compiled in DP) */
949 tmp2
= fr
->ener
[i
].esum
;
950 bOK
= bOK
&& gmx_fio_do_real(ef
->fio
, tmp2
);
952 fr
->ener
[i
].esum
= tmp2
;
954 if (file_version
== 1)
956 /* Old, unused real */
958 bOK
= bOK
&& gmx_fio_do_real(ef
->fio
, rdum
);
963 /* Here we can not check for file_version==1, since one could have
964 * continued an old format simulation with a new one with mdrun -append.
966 if (bRead
&& ef
->eo
.bOldFileOpen
)
968 /* Convert old full simulation sums to sums between energy frames */
969 convert_full_sums(&(ef
->eo
),fr
);
971 /* read the blocks */
972 for(b
=0; b
<fr
->nblock
; b
++)
974 /* now read the subblocks. */
975 int nsub
=fr
->block
[b
].nsub
; /* shortcut */
980 t_enxsubblock
*sub
=&(fr
->block
[b
].sub
[i
]); /* shortcut */
984 enxsubblock_alloc(sub
);
987 /* read/write data */
991 case xdr_datatype_float
:
992 bOK1
=gmx_fio_ndo_float(ef
->fio
, sub
->fval
, sub
->nr
);
994 case xdr_datatype_double
:
995 bOK1
=gmx_fio_ndo_double(ef
->fio
, sub
->dval
, sub
->nr
);
997 case xdr_datatype_int
:
998 bOK1
=gmx_fio_ndo_int(ef
->fio
, sub
->ival
, sub
->nr
);
1000 case xdr_datatype_large_int
:
1001 bOK1
=gmx_fio_ndo_gmx_large_int(ef
->fio
, sub
->lval
, sub
->nr
);
1003 case xdr_datatype_char
:
1004 bOK1
=gmx_fio_ndo_uchar(ef
->fio
, sub
->cval
, sub
->nr
);
1006 case xdr_datatype_string
:
1007 bOK1
=gmx_fio_ndo_string(ef
->fio
, sub
->sval
, sub
->nr
);
1010 gmx_incons("Reading unknown block data type: this file is corrupted or from the future");
1018 if( gmx_fio_flush(ef
->fio
) != 0)
1020 gmx_file("Cannot write energy file; maybe you are out of quota?");
1028 fprintf(stderr
,"\nLast energy frame read %d",
1030 fprintf(stderr
,"\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
1035 gmx_fatal(FARGS
,"could not write energies");
1043 static real
find_energy(const char *name
, int nre
, gmx_enxnm_t
*enm
,
1048 for(i
=0; i
<nre
; i
++)
1050 if (strcmp(enm
[i
].name
,name
) == 0)
1052 return fr
->ener
[i
].e
;
1056 gmx_fatal(FARGS
,"Could not find energy term named '%s'",name
);
1062 void get_enx_state(const char *fn
, real t
, gmx_groups_t
*groups
, t_inputrec
*ir
,
1065 /* Should match the names in mdebin.c */
1066 static const char *boxvel_nm
[] = {
1067 "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
1068 "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
1071 static const char *pcouplmu_nm
[] = {
1072 "Pcoupl-Mu-XX", "Pcoupl-Mu-YY", "Pcoupl-Mu-ZZ",
1073 "Pcoupl-Mu-YX", "Pcoupl-Mu-ZX", "Pcoupl-Mu-ZY"
1075 static const char *baro_nm
[] = {
1080 int ind0
[] = { XX
,YY
,ZZ
,YY
,ZZ
,ZZ
};
1081 int ind1
[] = { XX
,YY
,ZZ
,XX
,XX
,YY
};
1082 int nre
,nfr
,i
,j
,ni
,npcoupl
;
1085 gmx_enxnm_t
*enm
=NULL
;
1089 in
= open_enx(fn
,"r");
1090 do_enxnms(in
,&nre
,&enm
);
1093 while ((nfr
==0 || fr
->t
!= t
) && do_enx(in
,fr
)) {
1097 fprintf(stderr
,"\n");
1099 if (nfr
== 0 || fr
->t
!= t
)
1100 gmx_fatal(FARGS
,"Could not find frame with time %f in '%s'",t
,fn
);
1102 npcoupl
= TRICLINIC(ir
->compress
) ? 6 : 3;
1103 if (ir
->epc
== epcPARRINELLORAHMAN
) {
1104 clear_mat(state
->boxv
);
1105 for(i
=0; i
<npcoupl
; i
++) {
1106 state
->boxv
[ind0
[i
]][ind1
[i
]] =
1107 find_energy(boxvel_nm
[i
],nre
,enm
,fr
);
1109 fprintf(stderr
,"\nREAD %d BOX VELOCITIES FROM %s\n\n",npcoupl
,fn
);
1112 if (ir
->etc
== etcNOSEHOOVER
)
1118 for(i
=0; i
<state
->ngtc
; i
++) {
1119 ni
= groups
->grps
[egcTC
].nm_ind
[i
];
1120 bufi
= *(groups
->grpname
[ni
]);
1121 for(j
=0; (j
<state
->nhchainlength
); j
++)
1123 if (IR_NVT_TROTTER(ir
))
1125 sprintf(cns
,"-%d",j
);
1127 sprintf(buf
,"Xi%s-%s",cns
,bufi
);
1128 state
->nosehoover_xi
[i
] = find_energy(buf
,nre
,enm
,fr
);
1129 sprintf(buf
,"vXi%s-%s",cns
,bufi
);
1130 state
->nosehoover_vxi
[i
] = find_energy(buf
,nre
,enm
,fr
);
1134 fprintf(stderr
,"\nREAD %d NOSE-HOOVER Xi chains FROM %s\n\n",state
->ngtc
,fn
);
1136 if (IR_NPT_TROTTER(ir
))
1138 for(i
=0; i
<state
->nnhpres
; i
++) {
1139 bufi
= baro_nm
[0]; /* All barostat DOF's together for now */
1140 for(j
=0; (j
<state
->nhchainlength
); j
++)
1142 sprintf(buf
,"Xi-%d-%s",j
,bufi
);
1143 state
->nhpres_xi
[i
] = find_energy(buf
,nre
,enm
,fr
);
1144 sprintf(buf
,"vXi-%d-%s",j
,bufi
);
1145 state
->nhpres_vxi
[i
] = find_energy(buf
,nre
,enm
,fr
);
1148 fprintf(stderr
,"\nREAD %d NOSE-HOOVER BAROSTAT Xi chains FROM %s\n\n",state
->nnhpres
,fn
);
1152 free_enxnms(nre
,enm
);