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,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
47 #include "gromacs/fileio/gmxfio.h"
48 #include "gromacs/fileio/gmxfio_xdr.h"
49 #include "gromacs/fileio/xdrf.h"
50 #include "gromacs/math/functions.h"
51 #include "gromacs/math/vec.h"
52 #include "gromacs/mdtypes/inputrec.h"
53 #include "gromacs/mdtypes/md_enums.h"
54 #include "gromacs/mdtypes/state.h"
55 #include "gromacs/pbcutil/pbc.h"
56 #include "gromacs/topology/topology.h"
57 #include "gromacs/trajectory/energyframe.h"
58 #include "gromacs/utility/compare.h"
59 #include "gromacs/utility/fatalerror.h"
60 #include "gromacs/utility/futil.h"
61 #include "gromacs/utility/gmxassert.h"
62 #include "gromacs/utility/smalloc.h"
64 /* The source code in this file should be thread-safe.
65 Please keep it that way. */
67 /* This number should be increased whenever the file format changes! */
68 static const int enx_version = 5;
70 const char* enx_block_id_name[] = { "Averaged orientation restraints",
71 "Instantaneous orientation restraints",
72 "Orientation restraint order tensor(s)",
73 "Distance restraints",
80 /* Stuff for reading pre 4.1 energy files */
83 gmx_bool bOldFileOpen; /* Is this an open old file? */
84 gmx_bool bReadFirstStep; /* Did we read the first step? */
85 int first_step; /* First step in the energy file */
86 int step_prev; /* Previous step */
87 int nsum_prev; /* Previous step sum length */
88 t_energy* ener_prev; /* Previous energy sums */
99 static void enxsubblock_init(t_enxsubblock* sb)
103 sb->type = XdrDataType::Double;
105 sb->type = XdrDataType::Float;
121 static void enxsubblock_free(t_enxsubblock* sb)
157 for (i = 0; i < sb->sval_alloc; i++)
170 /* allocate the appropriate amount of memory for the given type and nr */
171 static void enxsubblock_alloc(t_enxsubblock* sb)
173 /* allocate the appropriate amount of memory */
176 case XdrDataType::Float:
177 if (sb->nr > sb->fval_alloc)
179 srenew(sb->fval, sb->nr);
180 sb->fval_alloc = sb->nr;
183 case XdrDataType::Double:
184 if (sb->nr > sb->dval_alloc)
186 srenew(sb->dval, sb->nr);
187 sb->dval_alloc = sb->nr;
190 case XdrDataType::Int:
191 if (sb->nr > sb->ival_alloc)
193 srenew(sb->ival, sb->nr);
194 sb->ival_alloc = sb->nr;
197 case XdrDataType::Int64:
198 if (sb->nr > sb->lval_alloc)
200 srenew(sb->lval, sb->nr);
201 sb->lval_alloc = sb->nr;
204 case XdrDataType::Char:
205 if (sb->nr > sb->cval_alloc)
207 srenew(sb->cval, sb->nr);
208 sb->cval_alloc = sb->nr;
211 case XdrDataType::String:
212 if (sb->nr > sb->sval_alloc)
216 srenew(sb->sval, sb->nr);
217 for (i = sb->sval_alloc; i < sb->nr; i++)
219 sb->sval[i] = nullptr;
221 sb->sval_alloc = sb->nr;
224 default: gmx_incons("Unknown block type: this file is corrupted or from the future");
228 static void enxblock_init(t_enxblock* eb)
236 static void enxblock_free(t_enxblock* eb)
238 if (eb->nsub_alloc > 0)
241 for (i = 0; i < eb->nsub_alloc; i++)
243 enxsubblock_free(&(eb->sub[i]));
251 void init_enxframe(t_enxframe* fr)
260 fr->nblock_alloc = 0;
265 void free_enxframe(t_enxframe* fr)
273 for (b = 0; b < fr->nblock_alloc; b++)
275 enxblock_free(&(fr->block[b]));
280 void add_blocks_enxframe(t_enxframe* fr, int n)
283 if (n > fr->nblock_alloc)
287 srenew(fr->block, n);
288 for (b = fr->nblock_alloc; b < fr->nblock; b++)
290 enxblock_init(&(fr->block[b]));
292 fr->nblock_alloc = n;
296 t_enxblock* find_block_id_enxframe(t_enxframe* ef, int id, t_enxblock* prev)
298 gmx_off_t starti = 0;
303 starti = (prev - ef->block) + 1;
305 for (i = starti; i < ef->nblock; i++)
307 if (ef->block[i].id == id)
309 return &(ef->block[i]);
315 void add_subblocks_enxblock(t_enxblock* eb, int n)
318 if (eb->nsub > eb->nsub_alloc)
323 for (b = eb->nsub_alloc; b < n; b++)
325 enxsubblock_init(&(eb->sub[b]));
331 static void enx_warning(const char* msg)
333 if (getenv("GMX_ENX_NO_FATAL") != nullptr)
335 gmx_warning("%s", msg);
342 "If you want to use the correct frames before the corrupted frame and avoid this "
343 "fatal error set the env.var. GMX_ENX_NO_FATAL");
347 static void edr_strings(XDR* xdr, gmx_bool bRead, int file_version, int n, gmx_enxnm_t** nms)
356 for (i = 0; i < n; i++)
372 if (!xdr_string(xdr, &(nm->name), STRLEN))
374 gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
376 if (file_version >= 2)
378 if (!xdr_string(xdr, &(nm->unit), STRLEN))
380 gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
385 nm->unit = gmx_strdup("kJ/mol");
390 void do_enxnms(ener_file_t ef, int* nre, gmx_enxnm_t** nms)
394 gmx_bool bRead = gmx_fio_getread(ef->fio);
397 xdr = gmx_fio_getxdr(ef->fio);
399 if (!xdr_int(xdr, &magic))
403 gmx_file("Cannot write energy names to file; maybe you are out of disk space?");
410 /* Assume this is an old edr format */
413 ef->eo.bOldFileOpen = TRUE;
414 ef->eo.bReadFirstStep = FALSE;
415 srenew(ef->eo.ener_prev, *nre);
419 ef->eo.bOldFileOpen = FALSE;
423 gmx_fatal(FARGS, "Energy names magic number mismatch, this is not a GROMACS edr file");
425 file_version = enx_version;
426 xdr_int(xdr, &file_version);
427 if (file_version > enx_version)
430 "reading tpx file (%s) version %d with version %d program",
431 gmx_fio_getname(ef->fio),
437 if (file_version != enx_version)
439 fprintf(stderr, "Note: enx file_version %d, software version %d\n", file_version, enx_version);
442 edr_strings(xdr, bRead, file_version, *nre, nms);
445 static gmx_bool do_eheader(ener_file_t ef,
449 gmx_bool* bWrongPrecision,
452 int magic = -7777777;
453 real first_real_to_check;
454 int b, zero = 0, dum = 0;
455 gmx_bool bRead = gmx_fio_getread(ef->fio);
459 XdrDataType xdrDataType = XdrDataType::Float;
461 XdrDataType xdrDataType = XdrDataType::Double;
466 *bWrongPrecision = FALSE;
470 /* The original energy frame started with a real,
471 * so we have to use a real for compatibility.
472 * This is VERY DIRTY code, since do_eheader can be called
473 * with the wrong precision set and then we could read r > -1e10,
474 * while actually the intention was r < -1e10.
475 * When nre_test >= 0, do_eheader should therefore terminate
476 * before the number of i/o calls starts depending on what has been read
477 * (which is the case for for instance the block sizes for variable
478 * number of blocks, where this number is read before).
480 first_real_to_check = -2e10;
481 if (!gmx_fio_do_real(ef->fio, first_real_to_check))
485 if (first_real_to_check > -1e10)
487 /* Assume we are reading an old format */
489 fr->t = first_real_to_check;
490 if (!gmx_fio_do_int(ef->fio, dum))
498 if (!gmx_fio_do_int(ef->fio, magic))
502 if (magic != -7777777)
504 enx_warning("Energy header magic number mismatch, this is not a GROMACS edr file");
508 *file_version = enx_version;
509 if (!gmx_fio_do_int(ef->fio, *file_version))
513 if (*bOK && *file_version > enx_version)
516 "reading tpx file (%s) version %d with version %d program",
517 gmx_fio_getname(ef->fio),
521 if (!gmx_fio_do_double(ef->fio, fr->t))
525 if (!gmx_fio_do_int64(ef->fio, fr->step))
529 if (!bRead && fr->nsum == 1)
531 /* Do not store sums of length 1,
532 * since this does not add information.
534 if (!gmx_fio_do_int(ef->fio, zero))
541 if (!gmx_fio_do_int(ef->fio, fr->nsum))
546 if (*file_version >= 3)
548 if (!gmx_fio_do_int64(ef->fio, fr->nsteps))
555 fr->nsteps = std::max(1, fr->nsum);
557 if (*file_version >= 5)
559 if (!gmx_fio_do_double(ef->fio, fr->dt))
569 if (!gmx_fio_do_int(ef->fio, fr->nre))
573 if (*file_version < 4)
575 if (!gmx_fio_do_int(ef->fio, ndisre))
582 /* now reserved for possible future use */
583 if (!gmx_fio_do_int(ef->fio, dum))
589 if (!gmx_fio_do_int(ef->fio, fr->nblock))
600 if (*file_version >= 4)
602 enx_warning("Distance restraint blocks in old style in new style file");
610 /* Frames could have nre=0, so we can not rely only on the fr->nre check */
611 if (bRead && nre_test >= 0
612 && ((fr->nre > 0 && fr->nre != nre_test) || fr->nre < 0 || ndisre < 0 || fr->nblock < 0))
616 *bWrongPrecision = TRUE;
621 /* we now know what these should be, or we've already bailed out because
622 of wrong precision */
623 if (*file_version == 1 && (fr->t < 0 || fr->t > 1e20 || fr->step < 0))
626 "edr file with negative step number or unreasonable time (and without version "
635 add_blocks_enxframe(fr, fr->nblock);
641 /* sub[0] is the instantaneous data, sub[1] is time averaged */
642 add_subblocks_enxblock(&(fr->block[0]), 2);
643 fr->block[0].id = enxDISRE;
644 fr->block[0].sub[0].nr = ndisre;
645 fr->block[0].sub[1].nr = ndisre;
646 fr->block[0].sub[0].type = xdrDataType;
647 fr->block[0].sub[1].type = xdrDataType;
651 /* read block header info */
652 for (b = startb; b < fr->nblock; b++)
654 if (*file_version < 4)
656 /* blocks in old version files always have 1 subblock that
657 consists of reals. */
662 add_subblocks_enxblock(&(fr->block[b]), 1);
666 if (fr->block[b].nsub != 1)
668 gmx_incons("Writing an old version .edr file with too many subblocks");
670 if (fr->block[b].sub[0].type != xdrDataType)
672 gmx_incons("Writing an old version .edr file the wrong subblock type");
675 nrint = fr->block[b].sub[0].nr;
677 if (!gmx_fio_do_int(ef->fio, nrint))
681 fr->block[b].id = b - startb;
682 fr->block[b].sub[0].nr = nrint;
683 fr->block[b].sub[0].type = xdrDataType;
688 /* in the new version files, the block header only contains
689 the ID and the number of subblocks */
690 int nsub = fr->block[b].nsub;
691 *bOK = *bOK && gmx_fio_do_int(ef->fio, fr->block[b].id);
692 *bOK = *bOK && gmx_fio_do_int(ef->fio, nsub);
694 fr->block[b].nsub = nsub;
697 add_subblocks_enxblock(&(fr->block[b]), nsub);
700 /* read/write type & size for each subblock */
701 for (i = 0; i < nsub; i++)
703 t_enxsubblock* sub = &(fr->block[b].sub[i]); /* shortcut */
704 int typenr = static_cast<int>(sub->type);
706 *bOK = *bOK && gmx_fio_do_int(ef->fio, typenr);
707 *bOK = *bOK && gmx_fio_do_int(ef->fio, sub->nr);
709 sub->type = static_cast<XdrDataType>(typenr);
713 if (!gmx_fio_do_int(ef->fio, fr->e_size))
718 /* now reserved for possible future use */
719 if (!gmx_fio_do_int(ef->fio, dum))
724 /* Do a dummy int to keep the format compatible with the old code */
725 if (!gmx_fio_do_int(ef->fio, dum))
730 if (*bOK && *file_version == 1 && nre_test < 0)
732 if (!ef->eo.bReadFirstStep)
734 ef->eo.bReadFirstStep = TRUE;
735 ef->eo.first_step = fr->step;
736 ef->eo.step_prev = fr->step;
737 ef->eo.nsum_prev = 0;
740 fr->nsum = fr->step - ef->eo.first_step + 1;
741 fr->nsteps = fr->step - ef->eo.step_prev;
748 void free_enxnms(int n, gmx_enxnm_t* nms)
752 for (i = 0; i < n; i++)
761 void close_enx(ener_file_t ef)
768 if (gmx_fio_close(ef->fio) != 0)
771 "Cannot close energy file; it might be corrupt, or maybe you are out of disk "
776 void done_ener_file(ener_file_t ef)
778 // Free the contents, then the pointer itself
783 /*!\brief Return TRUE if a file exists but is empty, otherwise FALSE.
785 * If the file exists but has length larger than zero, if it does not exist, or
786 * if there is a file system error, FALSE will be returned instead.
788 * \param fn File name to test
790 * \return TRUE if file could be open but is empty, otherwise FALSE.
792 static gmx_bool empty_file(const char* fn)
799 fp = gmx_fio_fopen(fn, "r");
800 ret = fread(&dum, sizeof(dum), 1, fp);
801 bEmpty = (feof(fp) != 0);
804 // bEmpty==TRUE but ret!=0 would likely be some strange I/O error, but at
805 // least it's not a normal empty file, so we return FALSE in that case.
806 return (bEmpty && ret == 0);
810 ener_file_t open_enx(const char* fn, const char* mode)
813 gmx_enxnm_t* nms = nullptr;
814 int file_version = -1;
816 gmx_bool bWrongPrecision, bOK = TRUE;
817 struct ener_file* ef;
823 ef->fio = gmx_fio_open(fn, mode);
824 gmx_fio_setprecision(ef->fio, FALSE);
825 do_enxnms(ef, &nre, &nms);
827 do_eheader(ef, &file_version, fr, nre, &bWrongPrecision, &bOK);
830 gmx_file("Cannot read energy file header. Corrupt file?");
833 /* Now check whether this file is in single precision */
835 && ((fr->e_size && (fr->nre == nre)
836 && (nre * 4 * static_cast<long int>(sizeof(float)) == fr->e_size))))
838 fprintf(stderr, "Opened %s as single precision energy file\n", fn);
839 free_enxnms(nre, nms);
843 gmx_fio_rewind(ef->fio);
844 gmx_fio_setprecision(ef->fio, TRUE);
845 do_enxnms(ef, &nre, &nms);
846 do_eheader(ef, &file_version, fr, nre, &bWrongPrecision, &bOK);
849 gmx_file("Cannot write energy file header; maybe you are out of disk space?");
852 if (((fr->e_size && (fr->nre == nre)
853 && (nre * 4 * static_cast<long int>(sizeof(double)) == fr->e_size))))
855 fprintf(stderr, "Opened %s as double precision energy file\n", fn);
861 gmx_fatal(FARGS, "File %s is empty", fn);
865 gmx_fatal(FARGS, "Energy file %s not recognized, maybe different CPU?", fn);
868 free_enxnms(nre, nms);
872 gmx_fio_rewind(ef->fio);
876 ef->fio = gmx_fio_open(fn, mode);
884 t_fileio* enx_file_pointer(const ener_file* ef)
889 static void convert_full_sums(ener_old_t* ener_old, t_enxframe* fr)
893 double esum_all, eav_all;
899 for (i = 0; i < fr->nre; i++)
901 if (fr->ener[i].e != 0)
905 if (fr->ener[i].esum != 0)
910 if (ne > 0 && ns == 0)
912 /* We do not have all energy sums */
917 /* Convert old full simulation sums to sums between energy frames */
918 nstep_all = fr->step - ener_old->first_step + 1;
919 if (fr->nsum > 1 && fr->nsum == nstep_all && ener_old->nsum_prev > 0)
921 /* Set the new sum length: the frame step difference */
922 fr->nsum = fr->step - ener_old->step_prev;
923 for (i = 0; i < fr->nre; i++)
925 esum_all = fr->ener[i].esum;
926 eav_all = fr->ener[i].eav;
927 fr->ener[i].esum = esum_all - ener_old->ener_prev[i].esum;
929 eav_all - ener_old->ener_prev[i].eav
930 - gmx::square(ener_old->ener_prev[i].esum / (nstep_all - fr->nsum) - esum_all / nstep_all)
931 * (nstep_all - fr->nsum) * nstep_all / static_cast<double>(fr->nsum);
932 ener_old->ener_prev[i].esum = esum_all;
933 ener_old->ener_prev[i].eav = eav_all;
935 ener_old->nsum_prev = nstep_all;
937 else if (fr->nsum > 0)
939 if (fr->nsum != nstep_all)
942 "\nWARNING: something is wrong with the energy sums, will not use exact "
944 ener_old->nsum_prev = 0;
948 ener_old->nsum_prev = nstep_all;
950 /* Copy all sums to ener_prev */
951 for (i = 0; i < fr->nre; i++)
953 ener_old->ener_prev[i].esum = fr->ener[i].esum;
954 ener_old->ener_prev[i].eav = fr->ener[i].eav;
958 ener_old->step_prev = fr->step;
961 gmx_bool do_enx(ener_file_t ef, t_enxframe* fr)
963 int file_version = -1;
965 gmx_bool bRead, bOK, bOK1, bSane;
966 real tmp1, tmp2, rdum;
970 bRead = gmx_fio_getread(ef->fio);
973 fr->e_size = fr->nre * sizeof(fr->ener[0].e) * 4;
974 /*d_size = fr->ndisre*(sizeof(real)*2);*/
977 if (!do_eheader(ef, &file_version, fr, -1, nullptr, &bOK))
981 fprintf(stderr, "\rLast energy frame read %d time %8.3f ", ef->framenr - 1, ef->frametime);
986 fprintf(stderr, "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n", ef->framenr, fr->t);
991 gmx_file("Cannot write energy file header; maybe you are out of disk space?");
997 if ((ef->framenr < 20 || ef->framenr % 10 == 0) && (ef->framenr < 200 || ef->framenr % 100 == 0)
998 && (ef->framenr < 2000 || ef->framenr % 1000 == 0))
1000 fprintf(stderr, "\rReading energy frame %6d time %8.3f ", ef->framenr, fr->t);
1003 ef->frametime = fr->t;
1005 /* Check sanity of this header */
1006 bSane = fr->nre > 0;
1007 for (b = 0; b < fr->nblock; b++)
1009 bSane = bSane || (fr->block[b].nsub > 0);
1011 if (!((fr->step >= 0) && bSane) && bRead)
1014 "\nWARNING: there may be something wrong with energy file %s\n",
1015 gmx_fio_getname(ef->fio));
1017 "Found: step=%" PRId64 ", nre=%d, nblock=%d, time=%g.\n",
1023 if (bRead && fr->nre > fr->e_alloc)
1025 srenew(fr->ener, fr->nre);
1026 for (i = fr->e_alloc; (i < fr->nre); i++)
1029 fr->ener[i].eav = 0;
1030 fr->ener[i].esum = 0;
1032 fr->e_alloc = fr->nre;
1035 for (i = 0; i < fr->nre; i++)
1037 bOK = bOK && gmx_fio_do_real(ef->fio, fr->ener[i].e);
1039 /* Do not store sums of length 1,
1040 * since this does not add information.
1042 if (file_version == 1 || (bRead && fr->nsum > 0) || fr->nsum > 1)
1044 tmp1 = fr->ener[i].eav;
1045 bOK = bOK && gmx_fio_do_real(ef->fio, tmp1);
1048 fr->ener[i].eav = tmp1;
1051 /* This is to save only in single precision (unless compiled in DP) */
1052 tmp2 = fr->ener[i].esum;
1053 bOK = bOK && gmx_fio_do_real(ef->fio, tmp2);
1056 fr->ener[i].esum = tmp2;
1059 if (file_version == 1)
1061 /* Old, unused real */
1063 bOK = bOK && gmx_fio_do_real(ef->fio, rdum);
1068 /* Here we can not check for file_version==1, since one could have
1069 * continued an old format simulation with a new one with mdrun -append.
1071 if (bRead && ef->eo.bOldFileOpen)
1073 /* Convert old full simulation sums to sums between energy frames */
1074 convert_full_sums(&(ef->eo), fr);
1076 /* read the blocks */
1077 for (b = 0; b < fr->nblock; b++)
1079 /* now read the subblocks. */
1080 int nsub = fr->block[b].nsub; /* shortcut */
1083 for (i = 0; i < nsub; i++)
1085 t_enxsubblock* sub = &(fr->block[b].sub[i]); /* shortcut */
1089 enxsubblock_alloc(sub);
1092 /* read/write data */
1095 case XdrDataType::Float:
1096 bOK1 = gmx_fio_ndo_float(ef->fio, sub->fval, sub->nr);
1098 case XdrDataType::Double:
1099 bOK1 = gmx_fio_ndo_double(ef->fio, sub->dval, sub->nr);
1101 case XdrDataType::Int: bOK1 = gmx_fio_ndo_int(ef->fio, sub->ival, sub->nr); break;
1102 case XdrDataType::Int64:
1103 bOK1 = gmx_fio_ndo_int64(ef->fio, sub->lval, sub->nr);
1105 case XdrDataType::Char:
1106 bOK1 = gmx_fio_ndo_uchar(ef->fio, sub->cval, sub->nr);
1108 case XdrDataType::String:
1109 bOK1 = gmx_fio_ndo_string(ef->fio, sub->sval, sub->nr);
1113 "Reading unknown block data type: this file is corrupted or from the "
1122 if (gmx_fio_flush(ef->fio) != 0)
1124 gmx_file("Cannot write energy file; maybe you are out of disk space?");
1132 fprintf(stderr, "\nLast energy frame read %d", ef->framenr - 1);
1133 fprintf(stderr, "\nWARNING: Incomplete energy frame: nr %d time %8.3f\n", ef->framenr, fr->t);
1137 gmx_fatal(FARGS, "could not write energies");
1145 static real find_energy(const char* name, int nre, gmx_enxnm_t* enm, t_enxframe* fr)
1149 for (i = 0; i < nre; i++)
1151 if (std::strcmp(enm[i].name, name) == 0)
1153 return fr->ener[i].e;
1158 "Could not find energy term named '%s'. Either the energy file is from a different "
1159 "run or this state variable is not stored in the energy file. In the latter case "
1160 "(and if you did not modify the T/P-coupling setup), you can read the state in mdrun "
1161 "instead, by passing in a checkpoint file.",
1166 void get_enx_state(const char* fn, real t, const SimulationGroups& groups, t_inputrec* ir, t_state* state)
1168 /* Should match the names in mdebin.c */
1169 static const char* boxvel_nm[] = { "Box-Vel-XX", "Box-Vel-YY", "Box-Vel-ZZ",
1170 "Box-Vel-YX", "Box-Vel-ZX", "Box-Vel-ZY" };
1172 static const char* baro_nm[] = { "Barostat" };
1175 int ind0[] = { XX, YY, ZZ, YY, ZZ, ZZ };
1176 int ind1[] = { XX, YY, ZZ, XX, XX, YY };
1177 int nre, nfr, i, j, ni, npcoupl;
1180 gmx_enxnm_t* enm = nullptr;
1184 in = open_enx(fn, "r");
1185 do_enxnms(in, &nre, &enm);
1188 while ((nfr == 0 || fr->t != t) && do_enx(in, fr))
1193 fprintf(stderr, "\n");
1195 if (nfr == 0 || fr->t != t)
1197 gmx_fatal(FARGS, "Could not find frame with time %f in '%s'", t, fn);
1200 npcoupl = TRICLINIC(ir->compress) ? 6 : 3;
1201 if (ir->epc == PressureCoupling::ParrinelloRahman)
1203 clear_mat(state->boxv);
1204 for (i = 0; i < npcoupl; i++)
1206 state->boxv[ind0[i]][ind1[i]] = find_energy(boxvel_nm[i], nre, enm, fr);
1208 fprintf(stderr, "\nREAD %d BOX VELOCITIES FROM %s\n\n", npcoupl, fn);
1211 if (ir->etc == TemperatureCoupling::NoseHoover)
1217 for (i = 0; i < state->ngtc; i++)
1219 ni = groups.groups[SimulationAtomGroupType::TemperatureCoupling][i];
1220 bufi = *(groups.groupNames[ni]);
1221 for (j = 0; (j < state->nhchainlength); j++)
1223 if (inputrecNvtTrotter(ir))
1225 sprintf(cns, "-%d", j);
1227 sprintf(buf, "Xi%s-%s", cns, bufi);
1228 state->nosehoover_xi[i] = find_energy(buf, nre, enm, fr);
1229 sprintf(buf, "vXi%s-%s", cns, bufi);
1230 state->nosehoover_vxi[i] = find_energy(buf, nre, enm, fr);
1233 fprintf(stderr, "\nREAD %d NOSE-HOOVER Xi chains FROM %s\n\n", state->ngtc, fn);
1235 if (inputrecNptTrotter(ir) || inputrecNphTrotter(ir))
1237 for (i = 0; i < state->nnhpres; i++)
1239 bufi = baro_nm[0]; /* All barostat DOF's together for now */
1240 for (j = 0; (j < state->nhchainlength); j++)
1242 sprintf(buf, "Xi-%d-%s", j, bufi);
1243 state->nhpres_xi[i] = find_energy(buf, nre, enm, fr);
1244 sprintf(buf, "vXi-%d-%s", j, bufi);
1245 state->nhpres_vxi[i] = find_energy(buf, nre, enm, fr);
1248 fprintf(stderr, "\nREAD %d NOSE-HOOVER BAROSTAT Xi chains FROM %s\n\n", state->nnhpres, fn);
1252 free_enxnms(nre, enm);
1257 static real ener_tensor_diag(int n,
1272 d1 = tensi[i] / DIM;
1273 d2 = tensi[i] - d1 * DIM;
1275 /* Find the diagonal elements d1 and d2 */
1276 len = std::strlen(enm1[ind1[i]].name);
1280 for (j = 0; j < n; j++)
1282 if (tensi[j] >= 0 && std::strlen(enm1[ind1[j]].name) == len
1283 && std::strncmp(enm1[ind1[i]].name, enm1[ind1[j]].name, len - 2) == 0
1284 && (tensi[j] == d1 * DIM + d1 || tensi[j] == d2 * DIM + d2))
1286 prod1 *= fabs(e1[ind1[j]].e);
1287 prod2 *= fabs(e2[ind2[j]].e);
1294 return 0.5 * (std::sqrt(prod1) + std::sqrt(prod2));
1302 static gmx_bool enernm_equal(const char* nm1, const char* nm2)
1306 len1 = std::strlen(nm1);
1307 len2 = std::strlen(nm2);
1309 /* Remove " (bar)" at the end of a name */
1310 if (len1 > 6 && std::strcmp(nm1 + len1 - 6, " (bar)") == 0)
1314 if (len2 > 6 && std::strcmp(nm2 + len2 - 6, " (bar)") == 0)
1319 return (len1 == len2 && gmx_strncasecmp(nm1, nm2, len1) == 0);
1322 static void cmp_energies(FILE* fp,
1336 int *tensi, len, d1, d2;
1337 real ftol_i, abstol_i;
1339 snew(tensi, maxener);
1340 /* Check for tensor elements ending on "-XX", "-XY", ... , "-ZZ" */
1341 for (i = 0; (i < maxener); i++)
1345 len = std::strlen(enm1[ii].name);
1346 if (len > 3 && enm1[ii].name[len - 3] == '-')
1348 d1 = enm1[ii].name[len - 2] - 'X';
1349 d2 = enm1[ii].name[len - 1] - 'X';
1350 if (d1 >= 0 && d1 < DIM && d2 >= 0 && d2 < DIM)
1352 tensi[i] = d1 * DIM + d2;
1357 for (i = 0; (i < maxener); i++)
1359 /* Check if this is an off-diagonal tensor element */
1360 if (tensi[i] >= 0 && tensi[i] != 0 && tensi[i] != 4 && tensi[i] != 8)
1362 /* Turn on the relative tolerance check (4 is maximum relative diff.) */
1364 /* Do the relative tolerance through an absolute tolerance times
1365 * the size of diagonal components of the tensor.
1367 abstol_i = ftol * ener_tensor_diag(nre, ind1, ind2, enm1, tensi, i, e1, e2);
1370 fprintf(debug, "tensor '%s' val %f diag %f\n", enm1[i].name, e1[i].e, abstol_i / ftol);
1374 /* We found a diagonal, we need to check with the minimum tolerance */
1375 abstol_i = std::min(abstol_i, abstol);
1379 /* We did not find a diagonal, ignore the relative tolerance check */
1388 if (!equal_real(e1[ind1[i]].e, e2[ind2[i]].e, ftol_i, abstol_i))
1391 "%-15s step %3d: %12g, step %3d: %12g\n",
1404 static void cmp_disres(t_enxframe *fr1, t_enxframe *fr2, real ftol, real abstol)
1407 char bav[64], bt[64], bs[22];
1409 cmp_int(stdout, "ndisre", -1, fr1->ndisre, fr2->ndisre);
1410 if ((fr1->ndisre == fr2->ndisre) && (fr1->ndisre > 0))
1412 sprintf(bav, "step %s: disre rav", gmx_step_str(fr1->step, bs));
1413 sprintf(bt, "step %s: disre rt", gmx_step_str(fr1->step, bs));
1414 for (i = 0; (i < fr1->ndisre); i++)
1416 cmp_real(stdout, bav, i, fr1->disre_rm3tav[i], fr2->disre_rm3tav[i], ftol, abstol);
1417 cmp_real(stdout, bt, i, fr1->disre_rt[i], fr2->disre_rt[i], ftol, abstol);
1423 static void cmp_eblocks(t_enxframe* fr1, t_enxframe* fr2, real ftol, real abstol)
1426 char buf[64], bs[22];
1428 cmp_int(stdout, "nblock", -1, fr1->nblock, fr2->nblock);
1429 if ((fr1->nblock == fr2->nblock) && (fr1->nblock > 0))
1431 for (j = 0; (j < fr1->nblock); j++)
1433 t_enxblock *b1, *b2; /* convenience vars */
1435 b1 = &(fr1->block[j]);
1436 b2 = &(fr2->block[j]);
1438 sprintf(buf, "step %s: block[%d]", gmx_step_str(fr1->step, bs), j);
1439 cmp_int(stdout, buf, -1, b1->nsub, b2->nsub);
1440 cmp_int(stdout, buf, -1, b1->id, b2->id);
1442 if ((b1->nsub == b2->nsub) && (b1->id == b2->id))
1444 for (i = 0; i < b1->nsub; i++)
1446 t_enxsubblock *s1, *s2;
1451 cmp_int(stdout, buf, -1, static_cast<int>(s1->type), static_cast<int>(s2->type));
1452 cmp_int64(stdout, buf, s1->nr, s2->nr);
1454 if ((s1->type == s2->type) && (s1->nr == s2->nr))
1458 case XdrDataType::Float:
1459 for (k = 0; k < s1->nr; k++)
1461 cmp_float(stdout, buf, i, s1->fval[k], s2->fval[k], ftol, abstol);
1464 case XdrDataType::Double:
1465 for (k = 0; k < s1->nr; k++)
1467 cmp_double(stdout, buf, i, s1->dval[k], s2->dval[k], ftol, abstol);
1470 case XdrDataType::Int:
1471 for (k = 0; k < s1->nr; k++)
1473 cmp_int(stdout, buf, i, s1->ival[k], s2->ival[k]);
1476 case XdrDataType::Int64:
1477 for (k = 0; k < s1->nr; k++)
1479 cmp_int64(stdout, buf, s1->lval[k], s2->lval[k]);
1482 case XdrDataType::Char:
1483 for (k = 0; k < s1->nr; k++)
1485 cmp_uc(stdout, buf, i, s1->cval[k], s2->cval[k]);
1488 case XdrDataType::String:
1489 for (k = 0; k < s1->nr; k++)
1491 cmp_str(stdout, buf, i, s1->sval[k], s2->sval[k]);
1494 default: gmx_incons("Unknown data type!!");
1503 void comp_enx(const char* fn1, const char* fn2, real ftol, real abstol, const char* lastener)
1505 int nre, nre1, nre2;
1506 ener_file_t in1, in2;
1507 int i, j, maxener, *ind1, *ind2, *have;
1508 gmx_enxnm_t *enm1 = nullptr, *enm2 = nullptr;
1509 t_enxframe * fr1, *fr2;
1512 fprintf(stdout, "comparing energy file %s and %s\n\n", fn1, fn2);
1514 in1 = open_enx(fn1, "r");
1515 in2 = open_enx(fn2, "r");
1516 do_enxnms(in1, &nre1, &enm1);
1517 do_enxnms(in2, &nre2, &enm2);
1520 fprintf(stdout, "There are %d and %d terms in the energy files\n\n", nre1, nre2);
1524 fprintf(stdout, "There are %d terms in the energy files\n\n", nre1);
1531 for (i = 0; i < nre1; i++)
1533 for (j = 0; j < nre2; j++)
1535 if (enernm_equal(enm1[i].name, enm2[j].name))
1544 if (nre == 0 || ind1[nre - 1] != i)
1546 cmp_str(stdout, "enm", i, enm1[i].name, "-");
1549 for (i = 0; i < nre2; i++)
1553 cmp_str(stdout, "enm", i, "-", enm2[i].name);
1558 for (i = 0; i < nre; i++)
1560 if ((lastener != nullptr) && (std::strstr(enm1[i].name, lastener) != nullptr))
1567 fprintf(stdout, "There are %d terms to compare in the energy files\n\n", maxener);
1569 for (i = 0; i < maxener; i++)
1571 cmp_str(stdout, "unit", i, enm1[ind1[i]].unit, enm2[ind2[i]].unit);
1578 b1 = do_enx(in1, fr1);
1579 b2 = do_enx(in2, fr2);
1582 fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn2, fn1);
1586 fprintf(stdout, "\nEnd of file on %s but not on %s\n", fn1, fn2);
1588 else if (!b1 && !b2)
1590 fprintf(stdout, "\nFiles read successfully\n");
1594 cmp_real(stdout, "t", -1, fr1->t, fr2->t, ftol, abstol);
1595 cmp_int(stdout, "step", -1, fr1->step, fr2->step);
1596 /* We don't want to print the nre mismatch for every frame */
1597 /* cmp_int(stdout,"nre",-1,fr1->nre,fr2->nre); */
1598 if ((fr1->nre >= nre) && (fr2->nre >= nre))
1601 stdout, fr1->step, fr1->step, fr1->ener, fr2->ener, enm1, ftol, abstol, nre, ind1, ind2, maxener);
1603 /*cmp_disres(fr1,fr2,ftol,abstol);*/
1604 cmp_eblocks(fr1, fr2, ftol, abstol);