2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015, 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.
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.
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.
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.
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.
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.
46 #include "gromacs/fileio/gmxfio.h"
47 #include "gromacs/fileio/gmxfio-xdr.h"
48 #include "gromacs/fileio/xdrf.h"
49 #include "gromacs/legacyheaders/macros.h"
50 #include "gromacs/math/vec.h"
51 #include "gromacs/pbcutil/pbc.h"
52 #include "gromacs/topology/topology.h"
53 #include "gromacs/utility/fatalerror.h"
54 #include "gromacs/utility/futil.h"
55 #include "gromacs/utility/gmxassert.h"
56 #include "gromacs/utility/smalloc.h"
58 /* The source code in this file should be thread-safe.
59 Please keep it that way. */
61 /* This number should be increased whenever the file format changes! */
62 static const int enx_version = 5;
64 const char *enx_block_id_name[] = {
65 "Averaged orientation restraints",
66 "Instantaneous orientation restraints",
67 "Orientation restraint order tensor(s)",
68 "Distance restraints",
75 /* Stuff for reading pre 4.1 energy files */
77 gmx_bool bOldFileOpen; /* Is this an open old file? */
78 gmx_bool bReadFirstStep; /* Did we read the first step? */
79 int first_step; /* First step in the energy file */
80 int step_prev; /* Previous step */
81 int nsum_prev; /* Previous step sum length */
82 t_energy *ener_prev; /* Previous energy sums */
93 static void enxsubblock_init(t_enxsubblock *sb)
97 sb->type = xdr_datatype_double;
99 sb->type = xdr_datatype_float;
115 static void enxsubblock_free(t_enxsubblock *sb)
151 for (i = 0; i < sb->sval_alloc; i++)
164 /* allocate the appropriate amount of memory for the given type and nr */
165 static void enxsubblock_alloc(t_enxsubblock *sb)
167 /* allocate the appropriate amount of memory */
170 case xdr_datatype_float:
171 if (sb->nr > sb->fval_alloc)
173 srenew(sb->fval, sb->nr);
174 sb->fval_alloc = sb->nr;
177 case xdr_datatype_double:
178 if (sb->nr > sb->dval_alloc)
180 srenew(sb->dval, sb->nr);
181 sb->dval_alloc = sb->nr;
184 case xdr_datatype_int:
185 if (sb->nr > sb->ival_alloc)
187 srenew(sb->ival, sb->nr);
188 sb->ival_alloc = sb->nr;
191 case xdr_datatype_int64:
192 if (sb->nr > sb->lval_alloc)
194 srenew(sb->lval, sb->nr);
195 sb->lval_alloc = sb->nr;
198 case xdr_datatype_char:
199 if (sb->nr > sb->cval_alloc)
201 srenew(sb->cval, sb->nr);
202 sb->cval_alloc = sb->nr;
205 case xdr_datatype_string:
206 if (sb->nr > sb->sval_alloc)
210 srenew(sb->sval, sb->nr);
211 for (i = sb->sval_alloc; i < sb->nr; i++)
215 sb->sval_alloc = sb->nr;
219 gmx_incons("Unknown block type: this file is corrupted or from the future");
223 static void enxblock_init(t_enxblock *eb)
231 static void enxblock_free(t_enxblock *eb)
233 if (eb->nsub_alloc > 0)
236 for (i = 0; i < eb->nsub_alloc; i++)
238 enxsubblock_free(&(eb->sub[i]));
246 void init_enxframe(t_enxframe *fr)
257 fr->nblock_alloc = 0;
262 void free_enxframe(t_enxframe *fr)
270 for (b = 0; b < fr->nblock_alloc; b++)
272 enxblock_free(&(fr->block[b]));
277 void add_blocks_enxframe(t_enxframe *fr, int n)
280 if (n > fr->nblock_alloc)
284 srenew(fr->block, n);
285 for (b = fr->nblock_alloc; b < fr->nblock; b++)
287 enxblock_init(&(fr->block[b]));
289 fr->nblock_alloc = n;
293 t_enxblock *find_block_id_enxframe(t_enxframe *ef, int id, t_enxblock *prev)
295 gmx_off_t starti = 0;
300 starti = (prev - ef->block) + 1;
302 for (i = starti; i < ef->nblock; i++)
304 if (ef->block[i].id == id)
306 return &(ef->block[i]);
312 void add_subblocks_enxblock(t_enxblock *eb, int n)
315 if (eb->nsub > eb->nsub_alloc)
320 for (b = eb->nsub_alloc; b < n; b++)
322 enxsubblock_init(&(eb->sub[b]));
328 static void enx_warning(const char *msg)
330 if (getenv("GMX_ENX_NO_FATAL") != NULL)
336 gmx_fatal(FARGS, "%s\n%s",
338 "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");
342 static void edr_strings(XDR *xdr, gmx_bool bRead, int file_version,
343 int n, gmx_enxnm_t **nms)
352 for (i = 0; i < n; i++)
368 if (!xdr_string(xdr, &(nm->name), STRLEN))
370 gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
372 if (file_version >= 2)
374 if (!xdr_string(xdr, &(nm->unit), STRLEN))
376 gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
381 nm->unit = gmx_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);
393 xdr = gmx_fio_getxdr(ef->fio);
395 if (!xdr_int(xdr, &magic))
399 gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
406 /* Assume this is an old edr format */
409 ef->eo.bOldFileOpen = TRUE;
410 ef->eo.bReadFirstStep = FALSE;
411 srenew(ef->eo.ener_prev, *nre);
415 ef->eo.bOldFileOpen = FALSE;
419 gmx_fatal(FARGS, "Energy names magic number mismatch, this is not a GROMACS edr file");
421 file_version = enx_version;
422 xdr_int(xdr, &file_version);
423 if (file_version > enx_version)
425 gmx_fatal(FARGS, "reading tpx file (%s) version %d with version %d program", gmx_fio_getname(ef->fio), file_version, enx_version);
429 if (file_version != enx_version)
431 fprintf(stderr, "Note: enx file_version %d, software version %d\n",
432 file_version, enx_version);
435 edr_strings(xdr, bRead, file_version, *nre, nms);
438 static gmx_bool do_eheader(ener_file_t ef, int *file_version, t_enxframe *fr,
439 int nre_test, gmx_bool *bWrongPrecision, gmx_bool *bOK)
441 int magic = -7777777;
442 real first_real_to_check;
443 int b, zero = 0, dum = 0;
444 gmx_bool bRead = gmx_fio_getread(ef->fio);
448 xdr_datatype dtreal = xdr_datatype_float;
450 xdr_datatype dtreal = xdr_datatype_double;
455 *bWrongPrecision = FALSE;
459 /* The original energy frame started with a real,
460 * so we have to use a real for compatibility.
461 * This is VERY DIRTY code, since do_eheader can be called
462 * with the wrong precision set and then we could read r > -1e10,
463 * while actually the intention was r < -1e10.
464 * When nre_test >= 0, do_eheader should therefore terminate
465 * before the number of i/o calls starts depending on what has been read
466 * (which is the case for for instance the block sizes for variable
467 * number of blocks, where this number is read before).
469 first_real_to_check = -2e10;
470 if (!gmx_fio_do_real(ef->fio, first_real_to_check))
474 if (first_real_to_check > -1e10)
476 /* Assume we are reading an old format */
478 fr->t = first_real_to_check;
479 if (!gmx_fio_do_int(ef->fio, dum))
487 if (!gmx_fio_do_int(ef->fio, magic))
491 if (magic != -7777777)
493 enx_warning("Energy header magic number mismatch, this is not a GROMACS edr file");
497 *file_version = enx_version;
498 if (!gmx_fio_do_int(ef->fio, *file_version))
502 if (*bOK && *file_version > enx_version)
504 gmx_fatal(FARGS, "reading tpx file (%s) version %d with version %d program", gmx_fio_getname(ef->fio), file_version, enx_version);
506 if (!gmx_fio_do_double(ef->fio, fr->t))
510 if (!gmx_fio_do_int64(ef->fio, fr->step))
514 if (!bRead && fr->nsum == 1)
516 /* Do not store sums of length 1,
517 * since this does not add information.
519 if (!gmx_fio_do_int(ef->fio, zero))
526 if (!gmx_fio_do_int(ef->fio, fr->nsum))
531 if (*file_version >= 3)
533 if (!gmx_fio_do_int64(ef->fio, fr->nsteps))
540 fr->nsteps = std::max(1, fr->nsum);
542 if (*file_version >= 5)
544 if (!gmx_fio_do_double(ef->fio, fr->dt))
554 if (!gmx_fio_do_int(ef->fio, fr->nre))
558 if (*file_version < 4)
560 if (!gmx_fio_do_int(ef->fio, ndisre))
567 /* now reserved for possible future use */
568 if (!gmx_fio_do_int(ef->fio, dum))
574 if (!gmx_fio_do_int(ef->fio, fr->nblock))
585 if (*file_version >= 4)
587 enx_warning("Distance restraint blocks in old style in new style file");
595 /* Frames could have nre=0, so we can not rely only on the fr->nre check */
596 if (bRead && nre_test >= 0 &&
597 ((fr->nre > 0 && fr->nre != nre_test) ||
598 fr->nre < 0 || ndisre < 0 || fr->nblock < 0))
602 *bWrongPrecision = TRUE;
607 /* we now know what these should be, or we've already bailed out because
608 of wrong precision */
609 if (*file_version == 1 && (fr->t < 0 || fr->t > 1e20 || fr->step < 0 ) )
611 enx_warning("edr file with negative step number or unreasonable time (and without version number).");
619 add_blocks_enxframe(fr, fr->nblock);
625 /* sub[0] is the instantaneous data, sub[1] is time averaged */
626 add_subblocks_enxblock(&(fr->block[0]), 2);
627 fr->block[0].id = enxDISRE;
628 fr->block[0].sub[0].nr = ndisre;
629 fr->block[0].sub[1].nr = ndisre;
630 fr->block[0].sub[0].type = dtreal;
631 fr->block[0].sub[1].type = dtreal;
635 /* read block header info */
636 for (b = startb; b < fr->nblock; b++)
638 if (*file_version < 4)
640 /* blocks in old version files always have 1 subblock that
641 consists of reals. */
646 add_subblocks_enxblock(&(fr->block[b]), 1);
650 if (fr->block[b].nsub != 1)
652 gmx_incons("Writing an old version .edr file with too many subblocks");
654 if (fr->block[b].sub[0].type != dtreal)
656 gmx_incons("Writing an old version .edr file the wrong subblock type");
659 nrint = fr->block[b].sub[0].nr;
661 if (!gmx_fio_do_int(ef->fio, nrint))
665 fr->block[b].id = b - startb;
666 fr->block[b].sub[0].nr = nrint;
667 fr->block[b].sub[0].type = dtreal;
672 /* in the new version files, the block header only contains
673 the ID and the number of subblocks */
674 int nsub = fr->block[b].nsub;
675 *bOK = *bOK && gmx_fio_do_int(ef->fio, fr->block[b].id);
676 *bOK = *bOK && gmx_fio_do_int(ef->fio, nsub);
678 fr->block[b].nsub = nsub;
681 add_subblocks_enxblock(&(fr->block[b]), nsub);
684 /* read/write type & size for each subblock */
685 for (i = 0; i < nsub; i++)
687 t_enxsubblock *sub = &(fr->block[b].sub[i]); /* shortcut */
688 int typenr = sub->type;
690 *bOK = *bOK && gmx_fio_do_int(ef->fio, typenr);
691 *bOK = *bOK && gmx_fio_do_int(ef->fio, sub->nr);
693 sub->type = (xdr_datatype)typenr;
697 if (!gmx_fio_do_int(ef->fio, fr->e_size))
702 /* now reserved for possible future use */
703 if (!gmx_fio_do_int(ef->fio, dum))
708 /* Do a dummy int to keep the format compatible with the old code */
709 if (!gmx_fio_do_int(ef->fio, dum))
714 if (*bOK && *file_version == 1 && nre_test < 0)
716 if (!ef->eo.bReadFirstStep)
718 ef->eo.bReadFirstStep = TRUE;
719 ef->eo.first_step = fr->step;
720 ef->eo.step_prev = fr->step;
721 ef->eo.nsum_prev = 0;
724 fr->nsum = fr->step - ef->eo.first_step + 1;
725 fr->nsteps = fr->step - ef->eo.step_prev;
732 void free_enxnms(int n, gmx_enxnm_t *nms)
736 for (i = 0; i < n; i++)
745 void close_enx(ener_file_t ef)
747 if (gmx_fio_close(ef->fio) != 0)
749 gmx_file("Cannot close energy file; it might be corrupt, or maybe you are out of disk space?");
753 /*!\brief Return TRUE if a file exists but is empty, otherwise FALSE.
755 * If the file exists but has length larger than zero, if it does not exist, or
756 * if there is a file system error, FALSE will be returned instead.
758 * \param fn File name to test
760 * \return TRUE if file could be open but is empty, otherwise FALSE.
763 empty_file(const char *fn)
770 fp = gmx_fio_fopen(fn, "r");
771 ret = fread(&dum, sizeof(dum), 1, fp);
775 // bEmpty==TRUE but ret!=0 would likely be some strange I/O error, but at
776 // least it's not a normal empty file, so we return FALSE in that case.
777 return (bEmpty && ret == 0);
781 ener_file_t open_enx(const char *fn, const char *mode)
784 gmx_enxnm_t *nms = NULL;
785 int file_version = -1;
787 gmx_bool bWrongPrecision, bOK = TRUE;
788 struct ener_file *ef;
794 ef->fio = gmx_fio_open(fn, mode);
795 gmx_fio_setprecision(ef->fio, FALSE);
796 do_enxnms(ef, &nre, &nms);
798 do_eheader(ef, &file_version, fr, nre, &bWrongPrecision, &bOK);
801 gmx_file("Cannot read energy file header. Corrupt file?");
804 /* Now check whether this file is in single precision */
805 if (!bWrongPrecision &&
806 ((fr->e_size && (fr->nre == nre) &&
807 (nre*4*(long int)sizeof(float) == fr->e_size)) ) )
809 fprintf(stderr, "Opened %s as single precision energy file\n", fn);
810 free_enxnms(nre, nms);
814 gmx_fio_rewind(ef->fio);
815 gmx_fio_setprecision(ef->fio, TRUE);
816 do_enxnms(ef, &nre, &nms);
817 do_eheader(ef, &file_version, fr, nre, &bWrongPrecision, &bOK);
820 gmx_file("Cannot write energy file header; maybe you are out of disk space?");
823 if (((fr->e_size && (fr->nre == nre) &&
824 (nre*4*(long int)sizeof(double) == fr->e_size)) ))
826 fprintf(stderr, "Opened %s as double precision energy file\n",
833 gmx_fatal(FARGS, "File %s is empty", fn);
837 gmx_fatal(FARGS, "Energy file %s not recognized, maybe different CPU?",
841 free_enxnms(nre, nms);
845 gmx_fio_rewind(ef->fio);
849 ef->fio = gmx_fio_open(fn, mode);
857 t_fileio *enx_file_pointer(const ener_file_t ef)
862 static void convert_full_sums(ener_old_t *ener_old, t_enxframe *fr)
866 double esum_all, eav_all;
872 for (i = 0; i < fr->nre; i++)
874 if (fr->ener[i].e != 0)
878 if (fr->ener[i].esum != 0)
883 if (ne > 0 && ns == 0)
885 /* We do not have all energy sums */
890 /* Convert old full simulation sums to sums between energy frames */
891 nstep_all = fr->step - ener_old->first_step + 1;
892 if (fr->nsum > 1 && fr->nsum == nstep_all && ener_old->nsum_prev > 0)
894 /* Set the new sum length: the frame step difference */
895 fr->nsum = fr->step - ener_old->step_prev;
896 for (i = 0; i < fr->nre; i++)
898 esum_all = fr->ener[i].esum;
899 eav_all = fr->ener[i].eav;
900 fr->ener[i].esum = esum_all - ener_old->ener_prev[i].esum;
901 fr->ener[i].eav = eav_all - ener_old->ener_prev[i].eav
902 - dsqr(ener_old->ener_prev[i].esum/(nstep_all - fr->nsum)
903 - esum_all/nstep_all)*
904 (nstep_all - fr->nsum)*nstep_all/(double)fr->nsum;
905 ener_old->ener_prev[i].esum = esum_all;
906 ener_old->ener_prev[i].eav = eav_all;
908 ener_old->nsum_prev = nstep_all;
910 else if (fr->nsum > 0)
912 if (fr->nsum != nstep_all)
914 fprintf(stderr, "\nWARNING: something is wrong with the energy sums, will not use exact averages\n");
915 ener_old->nsum_prev = 0;
919 ener_old->nsum_prev = nstep_all;
921 /* Copy all sums to ener_prev */
922 for (i = 0; i < fr->nre; i++)
924 ener_old->ener_prev[i].esum = fr->ener[i].esum;
925 ener_old->ener_prev[i].eav = fr->ener[i].eav;
929 ener_old->step_prev = fr->step;
932 gmx_bool do_enx(ener_file_t ef, t_enxframe *fr)
934 int file_version = -1;
936 gmx_bool bRead, bOK, bOK1, bSane;
937 real tmp1, tmp2, rdum;
941 bRead = gmx_fio_getread(ef->fio);
944 fr->e_size = fr->nre*sizeof(fr->ener[0].e)*4;
945 /*d_size = fr->ndisre*(sizeof(real)*2);*/
948 if (!do_eheader(ef, &file_version, fr, -1, NULL, &bOK))
952 fprintf(stderr, "\rLast energy frame read %d time %8.3f ",
953 ef->framenr-1, ef->frametime);
957 "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
963 gmx_file("Cannot write energy file header; maybe you are out of disk space?");
969 if ((ef->framenr < 20 || ef->framenr % 10 == 0) &&
970 (ef->framenr < 200 || ef->framenr % 100 == 0) &&
971 (ef->framenr < 2000 || ef->framenr % 1000 == 0))
973 fprintf(stderr, "\rReading energy frame %6d time %8.3f ",
977 ef->frametime = fr->t;
979 /* Check sanity of this header */
981 for (b = 0; b < fr->nblock; b++)
983 bSane = bSane || (fr->block[b].nsub > 0);
985 if (!((fr->step >= 0) && bSane))
987 fprintf(stderr, "\nWARNING: there may be something wrong with energy file %s\n",
988 gmx_fio_getname(ef->fio));
989 fprintf(stderr, "Found: step=%" GMX_PRId64 ", nre=%d, nblock=%d, time=%g.\n"
990 "Trying to skip frame expect a crash though\n",
991 fr->step, fr->nre, fr->nblock, fr->t);
993 if (bRead && fr->nre > fr->e_alloc)
995 srenew(fr->ener, fr->nre);
996 for (i = fr->e_alloc; (i < fr->nre); i++)
1000 fr->ener[i].esum = 0;
1002 fr->e_alloc = fr->nre;
1005 for (i = 0; i < fr->nre; i++)
1007 bOK = bOK && gmx_fio_do_real(ef->fio, fr->ener[i].e);
1009 /* Do not store sums of length 1,
1010 * since this does not add information.
1012 if (file_version == 1 ||
1013 (bRead && fr->nsum > 0) || fr->nsum > 1)
1015 tmp1 = fr->ener[i].eav;
1016 bOK = bOK && gmx_fio_do_real(ef->fio, tmp1);
1019 fr->ener[i].eav = tmp1;
1022 /* This is to save only in single precision (unless compiled in DP) */
1023 tmp2 = fr->ener[i].esum;
1024 bOK = bOK && gmx_fio_do_real(ef->fio, tmp2);
1027 fr->ener[i].esum = tmp2;
1030 if (file_version == 1)
1032 /* Old, unused real */
1034 bOK = bOK && gmx_fio_do_real(ef->fio, rdum);
1039 /* Here we can not check for file_version==1, since one could have
1040 * continued an old format simulation with a new one with mdrun -append.
1042 if (bRead && ef->eo.bOldFileOpen)
1044 /* Convert old full simulation sums to sums between energy frames */
1045 convert_full_sums(&(ef->eo), fr);
1047 /* read the blocks */
1048 for (b = 0; b < fr->nblock; b++)
1050 /* now read the subblocks. */
1051 int nsub = fr->block[b].nsub; /* shortcut */
1054 for (i = 0; i < nsub; i++)
1056 t_enxsubblock *sub = &(fr->block[b].sub[i]); /* shortcut */
1060 enxsubblock_alloc(sub);
1063 /* read/write data */
1066 case xdr_datatype_float:
1067 bOK1 = gmx_fio_ndo_float(ef->fio, sub->fval, sub->nr);
1069 case xdr_datatype_double:
1070 bOK1 = gmx_fio_ndo_double(ef->fio, sub->dval, sub->nr);
1072 case xdr_datatype_int:
1073 bOK1 = gmx_fio_ndo_int(ef->fio, sub->ival, sub->nr);
1075 case xdr_datatype_int64:
1076 bOK1 = gmx_fio_ndo_int64(ef->fio, sub->lval, sub->nr);
1078 case xdr_datatype_char:
1079 bOK1 = gmx_fio_ndo_uchar(ef->fio, sub->cval, sub->nr);
1081 case xdr_datatype_string:
1082 bOK1 = gmx_fio_ndo_string(ef->fio, sub->sval, sub->nr);
1085 gmx_incons("Reading unknown block data type: this file is corrupted or from the future");
1093 if (gmx_fio_flush(ef->fio) != 0)
1095 gmx_file("Cannot write energy file; maybe you are out of disk space?");
1103 fprintf(stderr, "\nLast energy frame read %d",
1105 fprintf(stderr, "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n",
1106 ef->framenr, fr->t);
1110 gmx_fatal(FARGS, "could not write energies");
1118 static real find_energy(const char *name, int nre, gmx_enxnm_t *enm,
1123 for (i = 0; i < nre; i++)
1125 if (std::strcmp(enm[i].name, name) == 0)
1127 return fr->ener[i].e;
1131 gmx_fatal(FARGS, "Could not find energy term named '%s'", name);
1137 void get_enx_state(const char *fn, real t, gmx_groups_t *groups, t_inputrec *ir,
1140 /* Should match the names in mdebin.c */
1141 static const char *boxvel_nm[] = {
1142 "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
1143 "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY"
1146 static const char *baro_nm[] = {
1151 int ind0[] = { XX, YY, ZZ, YY, ZZ, ZZ };
1152 int ind1[] = { XX, YY, ZZ, XX, XX, YY };
1153 int nre, nfr, i, j, ni, npcoupl;
1156 gmx_enxnm_t *enm = NULL;
1160 in = open_enx(fn, "r");
1161 do_enxnms(in, &nre, &enm);
1164 while ((nfr == 0 || fr->t != t) && do_enx(in, fr))
1169 fprintf(stderr, "\n");
1171 if (nfr == 0 || fr->t != t)
1173 gmx_fatal(FARGS, "Could not find frame with time %f in '%s'", t, fn);
1176 npcoupl = TRICLINIC(ir->compress) ? 6 : 3;
1177 if (ir->epc == epcPARRINELLORAHMAN)
1179 clear_mat(state->boxv);
1180 for (i = 0; i < npcoupl; i++)
1182 state->boxv[ind0[i]][ind1[i]] =
1183 find_energy(boxvel_nm[i], nre, enm, fr);
1185 fprintf(stderr, "\nREAD %d BOX VELOCITIES FROM %s\n\n", npcoupl, fn);
1188 if (ir->etc == etcNOSEHOOVER)
1194 for (i = 0; i < state->ngtc; i++)
1196 ni = groups->grps[egcTC].nm_ind[i];
1197 bufi = *(groups->grpname[ni]);
1198 for (j = 0; (j < state->nhchainlength); j++)
1200 if (IR_NVT_TROTTER(ir))
1202 sprintf(cns, "-%d", j);
1204 sprintf(buf, "Xi%s-%s", cns, bufi);
1205 state->nosehoover_xi[i] = find_energy(buf, nre, enm, fr);
1206 sprintf(buf, "vXi%s-%s", cns, bufi);
1207 state->nosehoover_vxi[i] = find_energy(buf, nre, enm, fr);
1211 fprintf(stderr, "\nREAD %d NOSE-HOOVER Xi chains FROM %s\n\n", state->ngtc, fn);
1213 if (IR_NPT_TROTTER(ir) || IR_NPH_TROTTER(ir))
1215 for (i = 0; i < state->nnhpres; i++)
1217 bufi = baro_nm[0]; /* All barostat DOF's together for now */
1218 for (j = 0; (j < state->nhchainlength); j++)
1220 sprintf(buf, "Xi-%d-%s", j, bufi);
1221 state->nhpres_xi[i] = find_energy(buf, nre, enm, fr);
1222 sprintf(buf, "vXi-%d-%s", j, bufi);
1223 state->nhpres_vxi[i] = find_energy(buf, nre, enm, fr);
1226 fprintf(stderr, "\nREAD %d NOSE-HOOVER BAROSTAT Xi chains FROM %s\n\n", state->nnhpres, fn);
1230 free_enxnms(nre, enm);