From e1f1aabf9cbf72f368e0be3f9647698dd8ebf1b3 Mon Sep 17 00:00:00 2001 From: Mark Abraham Date: Thu, 17 Jan 2019 21:01:02 +0100 Subject: [PATCH] Remove copy_energy from energy-writing steps Used range-based for to avoid yet another copy of energy data between intermediate buffers before being written to the .edr file. Renamed some variable from index to entryIndex for clarity. This change is covered by the tests on t_mdebin functionality. Change-Id: Iecaf28c8af48a74abaf2805feb353c99ab613b10 --- src/gromacs/mdlib/ebin.cpp | 76 +++++++++++++++++++++++++++++++++++++------- src/gromacs/mdlib/ebin.h | 22 +++++++++---- src/gromacs/mdlib/mdebin.cpp | 21 +----------- 3 files changed, 82 insertions(+), 37 deletions(-) diff --git a/src/gromacs/mdlib/ebin.cpp b/src/gromacs/mdlib/ebin.cpp index 50d24b0d9c..3c32852e90 100644 --- a/src/gromacs/mdlib/ebin.cpp +++ b/src/gromacs/mdlib/ebin.cpp @@ -42,14 +42,18 @@ #include #include +#include + #include "gromacs/math/units.h" #include "gromacs/math/utilities.h" #include "gromacs/math/vec.h" #include "gromacs/topology/ifunc.h" #include "gromacs/trajectory/energyframe.h" +#include "gromacs/utility/basedefinitions.h" #include "gromacs/utility/cstringutil.h" #include "gromacs/utility/fatalerror.h" #include "gromacs/utility/smalloc.h" +#include "gromacs/utility/stringutil.h" t_ebin *mk_ebin() { @@ -129,19 +133,19 @@ int get_ebin_space(t_ebin *eb, int nener, const char *const enm[], const char *u return index; } -void add_ebin(t_ebin *eb, int index, int nener, const real ener[], gmx_bool bSum) +void add_ebin(t_ebin *eb, int entryIndex, int nener, const real ener[], gmx_bool bSum) { int i, m; double e, invmm, diff; t_energy *eg, *egs; - if ((index+nener > eb->nener) || (index < 0)) + if ((entryIndex+nener > eb->nener) || (entryIndex < 0)) { - gmx_fatal(FARGS, "%s-%d: Energies out of range: index=%d nener=%d maxener=%d", - __FILE__, __LINE__, index, nener, eb->nener); + gmx_fatal(FARGS, "%s-%d: Energies out of range: entryIndex=%d nener=%d maxener=%d", + __FILE__, __LINE__, entryIndex, nener, eb->nener); } - eg = &(eb->e[index]); + eg = &(eb->e[entryIndex]); for (i = 0; (i < nener); i++) { @@ -150,7 +154,7 @@ void add_ebin(t_ebin *eb, int index, int nener, const real ener[], gmx_bool bSum if (bSum) { - egs = &(eb->e_sim[index]); + egs = &(eb->e_sim[entryIndex]); m = eb->nsum; @@ -182,6 +186,56 @@ void add_ebin(t_ebin *eb, int index, int nener, const real ener[], gmx_bool bSum } } +// TODO It would be faster if this function was templated on both bSum +// and whether eb->nsum was zero, to lift the branches out of the loop +// over all possible energy terms, but that is true for a lot of the +// ebin and mdebin functionality, so we should do it all or nothing. +void add_ebin_indexed(t_ebin *eb, int entryIndex, gmx::ArrayRef shouldUse, + gmx::ArrayRef ener, gmx_bool bSum) +{ + + GMX_ASSERT(shouldUse.size() == ener.size(), "View sizes must match"); + GMX_ASSERT(entryIndex+std::count(shouldUse.begin(), shouldUse.end(), true) <= eb->nener, + gmx::formatString("Energies out of range: entryIndex=%d nener=%zu maxener=%d", + entryIndex, + std::count(shouldUse.begin(), shouldUse.end(), true), + eb->nener).c_str()); + GMX_ASSERT(entryIndex >= 0, "Must have non-negative entry"); + + const int m = eb->nsum; + const double invmm = (m == 0) ? 0 : (1.0/m)/(m+1.0); + t_energy *energyEntry = &(eb->e[entryIndex]); + t_energy *simEnergyEntry = &(eb->e_sim[entryIndex]); + auto shouldUseIter = shouldUse.begin(); + for (const auto &theEnergy : ener) + { + if (*shouldUseIter) + { + energyEntry->e = theEnergy; + if (bSum) + { + if (m == 0) + { + energyEntry->eav = 0; + energyEntry->esum = theEnergy; + simEnergyEntry->esum += theEnergy; + } + else + { + /* first update sigma, then sum */ + double diff = energyEntry->esum - m*theEnergy; + energyEntry->eav += diff*diff*invmm; + energyEntry->esum += theEnergy; + simEnergyEntry->esum += theEnergy; + } + ++simEnergyEntry; + } + ++energyEntry; + } + ++shouldUseIter; + } +} + void ebin_increase_count(t_ebin *eb, gmx_bool bSum) { eb->nsteps++; @@ -201,7 +255,7 @@ void reset_ebin_sums(t_ebin *eb) /* The actual sums are cleared when the next frame is stored */ } -void pr_ebin(FILE *fp, t_ebin *eb, int index, int nener, int nperline, +void pr_ebin(FILE *fp, t_ebin *eb, int entryIndex, int nener, int nperline, int prmode, gmx_bool bPrHead) { int i, j, i0; @@ -210,11 +264,11 @@ void pr_ebin(FILE *fp, t_ebin *eb, int index, int nener, int nperline, rc = 0; - if (index < 0 || index > eb->nener) + if (entryIndex < 0 || entryIndex > eb->nener) { - gmx_fatal(FARGS, "Invalid index in pr_ebin: %d", index); + gmx_fatal(FARGS, "Invalid entryIndex in pr_ebin: %d", entryIndex); } - int start = index; + int start = entryIndex; if (nener > eb->nener) { gmx_fatal(FARGS, "Invalid nener in pr_ebin: %d", nener); @@ -222,7 +276,7 @@ void pr_ebin(FILE *fp, t_ebin *eb, int index, int nener, int nperline, int end = eb->nener; if (nener != -1) { - end = index + nener; + end = entryIndex + nener; } for (i = start; (i < end) && rc >= 0; ) { diff --git a/src/gromacs/mdlib/ebin.h b/src/gromacs/mdlib/ebin.h index 3a52717748..e843438bf9 100644 --- a/src/gromacs/mdlib/ebin.h +++ b/src/gromacs/mdlib/ebin.h @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2018, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2018,2019, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -41,6 +41,7 @@ #include "gromacs/fileio/enxio.h" #include "gromacs/trajectory/energyframe.h" +#include "gromacs/utility/arrayref.h" #include "gromacs/utility/basedefinitions.h" /* This is a running averaging structure ('energy bin') for use during mdrun. */ @@ -70,17 +71,26 @@ int get_ebin_space(t_ebin *eb, int nener, const char *const enm[], const char *u /* Create space in the energy bin and register names. * The enm array must be static, because the contents are not copied, * but only the pointers. - * Function returns an index number that must be used in subsequent + * Function returns an entryIndex number that must be used in subsequent * calls to add_ebin. */ -void add_ebin(t_ebin *eb, int index, int nener, const real ener[], gmx_bool bSum); +void add_ebin(t_ebin *eb, int entryIndex, int nener, const real ener[], gmx_bool bSum); /* Add nener reals (eg. energies, box-lengths, pressures) to the - * energy bin at position index. + * energy bin at position entryIndex. * If bSum is TRUE then the reals are also added to the sum * and sum of squares. */ +/*! \brief Add values from \c ener to \c eb if the matching entry in + * shouldUse is true. + * + * Caller must ensure that \c shouldUse and \c ener to have the same + * size, and that \c eb has enough room for the number of true + * entries in \c shouldUse. */ +void add_ebin_indexed(t_ebin *eb, int entryIndex, gmx::ArrayRef shouldUse, + gmx::ArrayRef ener, gmx_bool bSum); + void ebin_increase_count(t_ebin *eb, gmx_bool bSum); /* Increase the counters for the sums. * This routine should be called AFTER all add_ebin calls for this step. @@ -96,12 +106,12 @@ void reset_ebin_sums(t_ebin *eb); * used only when eprAVER or eprRMS is set. If bPrHead than the * header is printed. * - * \c index and \c nener must be in [0, eb->nener), except that \c + * \c entryIndex and \c nener must be in [0, eb->nener), except that \c * nener -1 is interpreted as eb->nener. * * \todo Callers should be refactored pass eb->nener, rather than * us implement and rely on this special behaviour of -1. */ -void pr_ebin(FILE *fp, t_ebin *eb, int index, int nener, int nperline, +void pr_ebin(FILE *fp, t_ebin *eb, int entryIndex, int nener, int nperline, int prmode, gmx_bool bPrHead); #endif diff --git a/src/gromacs/mdlib/mdebin.cpp b/src/gromacs/mdlib/mdebin.cpp index 8fd29e2fa5..4070597090 100644 --- a/src/gromacs/mdlib/mdebin.cpp +++ b/src/gromacs/mdlib/mdebin.cpp @@ -890,23 +890,6 @@ extern FILE *open_dhdl(const char *filename, const t_inputrec *ir, return fp; } -static void copy_energy(t_mdebin *md, const real e[], real ecpy[]) -{ - int i, j; - - for (i = j = 0; (i < F_NRE); i++) - { - if (md->bEner[i]) - { - ecpy[j++] = e[i]; - } - } - if (j != md->f_nre) - { - gmx_incons("Number of energy terms wrong"); - } -} - void upd_mdebin(t_mdebin *md, gmx_bool bDoDHDL, gmx_bool bSum, @@ -929,7 +912,6 @@ void upd_mdebin(t_mdebin *md, real crmsd[2], tmp6[6]; real bs[NTRICLBOXS], vol, dens, pv, enthalpy; real eee[egNR]; - real ecopy[F_NRE]; double store_dhdl[efptNR]; real store_energy = 0; real tmp; @@ -938,8 +920,7 @@ void upd_mdebin(t_mdebin *md, * as an argument. This is because we sometimes need to write the box from * the last timestep to match the trajectory frames. */ - copy_energy(md, enerd->term, ecopy); - add_ebin(md->ebin, md->ie, md->f_nre, ecopy, bSum); + add_ebin_indexed(md->ebin, md->ie, gmx::ArrayRef(md->bEner), enerd->term, bSum); if (md->nCrmsd) { crmsd[0] = constr->rmsd(); -- 2.11.4.GIT