#include <cmath>
#include <cstring>
+#include <algorithm>
+
#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()
{
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++)
{
if (bSum)
{
- egs = &(eb->e_sim[index]);
+ egs = &(eb->e_sim[entryIndex]);
m = eb->nsum;
}
}
+// 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<bool> shouldUse,
+ gmx::ArrayRef<const real> 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++;
/* 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;
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);
int end = eb->nener;
if (nener != -1)
{
- end = index + nener;
+ end = entryIndex + nener;
}
for (i = start; (i < end) && rc >= 0; )
{
*
* 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.
#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. */
/* 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<bool> shouldUse,
+ gmx::ArrayRef<const real> 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.
* 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
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,
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;
* 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<gmx_bool>(md->bEner), enerd->term, bSum);
if (md->nCrmsd)
{
crmsd[0] = constr->rmsd();
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff