Write pull data to checkpoints to keep a proper average.
Add option whether to write instantaneous or average pull
coordinates and forces, respectively.
Fixes #2626
Change-Id: Id0268e00486028f497463dc13ae8d65a6c5df325
With this option the PBC reference atom is only used at initialization.
This can be of use when using large pull groups or groups with potentially
large relative movement of atoms.
+
+Enable output of average pull forces and positions
+"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""
+
+Normally the pull module writes instantaneous output of positions and forces, however
+now it can write the average of these values over the period since the last output.
+This works correctly even if a checkpoint restart intervened. This is enabled using the
+new options ``pull-fout-average`` and ``pull-xout-average``.
be located centrally in the group. Using the COM from the
previous step can be useful if one or more pull groups are large.
+.. mdp:: pull-xout-average
+
+ .. mdp-value:: no
+
+ Write the instantaneous coordinates for all the pulled groups.
+
+ .. mdp-value:: yes
+
+ Write the average coordinates (since last output) for all the
+ pulled groups. N.b., some analysis tools might expect instantaneous
+ pull output.
+
+.. mdp:: pull-fout-average
+
+ .. mdp-value:: no
+
+ Write the instantaneous force for all the pulled groups.
+
+ .. mdp-value:: yes
+
+ Write the average force (since last output) for all the
+ pulled groups. N.b., some analysis tools might expect instantaneous
+ pull output.
+
.. mdp:: pull-ngroups
(1)
#include "gromacs/mdtypes/inputrec.h"
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/observableshistory.h"
+#include "gromacs/mdtypes/pullhistory.h"
#include "gromacs/mdtypes/state.h"
#include "gromacs/mdtypes/swaphistory.h"
#include "gromacs/trajectory/trajectoryframe.h"
cptv_Unknown = 17, /**< Version before numbering scheme */
cptv_RemoveBuildMachineInformation, /**< remove functionality that makes mdrun builds non-reproducible */
cptv_ComPrevStepAsPullGroupReference, /**< Allow using COM of previous step as pull group PBC reference */
+ cptv_PullAverage, /**< Added possibility to output average pull force and position */
cptv_Count /**< the total number of cptv versions */
};
eenhNR
};
+enum {
+ epullhPULL_NUMCOORDINATES, epullhPULL_NUMGROUPS, epullhPULL_NUMVALUESINXSUM,
+ epullhPULL_NUMVALUESINFSUM, epullhNR
+};
+
+enum {
+ epullcoordh_VALUE_REF_SUM, epullcoordh_VALUE_SUM, epullcoordh_DR01_SUM,
+ epullcoordh_DR23_SUM, epullcoordh_DR45_SUM, epullcoordh_FSCAL_SUM,
+ epullcoordh_DYNAX_SUM, epullcoordh_NR
+};
+
+enum {
+ epullgrouph_X_SUM, epullgrouph_NR
+};
+
static const char *eenh_names[eenhNR] =
{
"energy_n", "energy_aver", "energy_sum", "energy_nsum",
"energy_delta_h_start_lambda"
};
+static const char *ePullhNames[epullhNR] =
+{
+ "pullhistory_numcoordinates", "pullhistory_numgroups", "pullhistory_numvaluesinxsum",
+ "pullhistory_numvaluesinfsum"
+};
+
/* free energy history variables -- need to be preserved over checkpoint */
enum {
edfhBEQUIL, edfhNATLAMBDA, edfhWLHISTO, edfhWLDELTA, edfhSUMWEIGHTS, edfhSUMDG, edfhSUMMINVAR, edfhSUMVAR,
energyHistory, //!< Energy observable statistics
freeEnergyHistory, //!< Free-energy state and observable statistics
accWeightHistogram, //!< Accelerated weight histogram method state
- pullState //!< COM of previous step.
+ pullState, //!< COM of previous step.
+ pullHistory //!< Pull history statistics (sums since last written output)
};
//! \brief Return the name of a checkpoint entry based on part and part entry
case StatePart::freeEnergyHistory: return edfh_names[ecpt];
case StatePart::accWeightHistogram: return eawhh_names[ecpt];
case StatePart::pullState: return epull_prev_step_com_names[ecpt];
+ case StatePart::pullHistory: return ePullhNames[ecpt];
}
return nullptr;
return 0;
}
-//! \brief Read/Write a std::vector.
+//! \brief Read/Write a std::vector, on read checks the number of elements matches \p numElements, if specified.
template <typename T>
static int doVector(XDR *xd, StatePart part, int ecpt, int sflags,
- std::vector<T> *vector, FILE *list)
+ std::vector<T> *vector, FILE *list, int numElements = -1)
{
- return doVectorLow<T>(xd, part, ecpt, sflags, -1, nullptr, nullptr, vector, list, CptElementType::real);
+ return doVectorLow<T>(xd, part, ecpt, sflags, numElements, nullptr, nullptr, vector, list, CptElementType::real);
}
//! \brief Read/Write an ArrayRef<real>.
int flags_state;
int flags_eks;
int flags_enh;
+ int flagsPullHistory;
int flags_dfh;
int flags_awhh;
int nED;
{
contents->flags_awhh = 0;
}
+
+ if (contents->file_version >= 18)
+ {
+ do_cpt_int_err(xd, "pull history flags", &contents->flagsPullHistory, list);
+ }
+ else
+ {
+ contents->flagsPullHistory = 0;
+ }
}
static int do_cpt_footer(XDR *xd, int file_version)
return ret;
}
+static int doCptPullCoordHist(XDR *xd, PullCoordinateHistory *pullCoordHist,
+ const StatePart part, FILE *list)
+{
+ int ret = 0;
+ int flags = 0;
+
+ flags |= ((1<<epullcoordh_VALUE_REF_SUM) | (1<<epullcoordh_VALUE_SUM) | (1<<epullcoordh_DR01_SUM) |
+ (1<<epullcoordh_DR23_SUM) | (1<<epullcoordh_DR45_SUM) | (1<<epullcoordh_FSCAL_SUM));
+
+ for (int i = 0; i < epullcoordh_NR && ret == 0; i++)
+ {
+ switch (i)
+ {
+ case epullcoordh_VALUE_REF_SUM: ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->valueRef), list); break;
+ case epullcoordh_VALUE_SUM: ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->value), list); break;
+ case epullcoordh_DR01_SUM:
+ for (int j = 0; j < DIM && ret == 0; j++)
+ {
+ ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dr01[j]), list);
+ }
+ break;
+ case epullcoordh_DR23_SUM:
+ for (int j = 0; j < DIM && ret == 0; j++)
+ {
+ ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dr23[j]), list);
+ }
+ break;
+ case epullcoordh_DR45_SUM:
+ for (int j = 0; j < DIM && ret == 0; j++)
+ {
+ ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dr45[j]), list);
+ }
+ break;
+ case epullcoordh_FSCAL_SUM: ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->scalarForce), list); break;
+ case epullcoordh_DYNAX_SUM:
+ for (int j = 0; j < DIM && ret == 0; j++)
+ {
+ ret = do_cpte_double(xd, part, i, flags, &(pullCoordHist->dynaX[j]), list);
+ }
+ break;
+ }
+ }
+
+ return ret;
+}
+
+static int doCptPullGroupHist(XDR *xd, PullGroupHistory *pullGroupHist,
+ const StatePart part, FILE *list)
+{
+ int ret = 0;
+ int flags = 0;
+
+ flags |= ((1<<epullgrouph_X_SUM));
+
+ for (int i = 0; i < epullgrouph_NR; i++)
+ {
+ switch (i)
+ {
+ case epullgrouph_X_SUM:
+ for (int j = 0; j < DIM && ret == 0; j++)
+ {
+ ret = do_cpte_double(xd, part, i, flags, &(pullGroupHist->x[j]), list);
+ }
+ break;
+ }
+ }
+
+ return ret;
+}
+
+
+static int doCptPullHist(XDR *xd, gmx_bool bRead,
+ int fflags, PullHistory *pullHist,
+ const StatePart part,
+ FILE *list)
+{
+ int ret = 0;
+ int pullHistoryNumCoordinates = 0;
+ int pullHistoryNumGroups = 0;
+
+ /* Retain the number of terms in the sum and the number of coordinates (used for writing
+ * average pull forces and coordinates) in the pull history, in temporary variables,
+ * in case they cannot be read from the checkpoint, in order to have backward compatibility */
+ if (bRead)
+ {
+ pullHist->numValuesInXSum = 0;
+ pullHist->numValuesInFSum = 0;
+ }
+ else if (pullHist != nullptr)
+ {
+ pullHistoryNumCoordinates = pullHist->pullCoordinateSums.size();
+ pullHistoryNumGroups = pullHist->pullGroupSums.size();
+ }
+ else
+ {
+ GMX_RELEASE_ASSERT(fflags == 0, "Without pull history, all flags should be off");
+ }
+
+ for (int i = 0; (i < epullhNR && ret == 0); i++)
+ {
+ if (fflags & (1<<i))
+ {
+ switch (i)
+ {
+ case epullhPULL_NUMCOORDINATES: ret = do_cpte_int(xd, part, i, fflags, &pullHistoryNumCoordinates, list); break;
+ case epullhPULL_NUMGROUPS: do_cpt_int_err(xd, eenh_names[i], &pullHistoryNumGroups, list); break;
+ case epullhPULL_NUMVALUESINXSUM: do_cpt_int_err(xd, eenh_names[i], &pullHist->numValuesInXSum, list); break;
+ case epullhPULL_NUMVALUESINFSUM: do_cpt_int_err(xd, eenh_names[i], &pullHist->numValuesInFSum, list); break;
+ default:
+ gmx_fatal(FARGS, "Unknown pull history entry %d\n"
+ "You are probably reading a new checkpoint file with old code", i);
+ }
+ }
+ }
+ if (bRead)
+ {
+ pullHist->pullCoordinateSums.resize(pullHistoryNumCoordinates);
+ pullHist->pullGroupSums.resize(pullHistoryNumGroups);
+ }
+ if (pullHist->numValuesInXSum > 0 || pullHist->numValuesInFSum > 0)
+ {
+ for (size_t i = 0; i < pullHist->pullCoordinateSums.size() && ret == 0; i++)
+ {
+ ret = doCptPullCoordHist(xd, &(pullHist->pullCoordinateSums[i]), part, list);
+ }
+ for (size_t i = 0; i < pullHist->pullGroupSums.size() && ret == 0; i++)
+ {
+ ret = doCptPullGroupHist(xd, &(pullHist->pullGroupSums[i]), part, list);
+ }
+ }
+
+ return ret;
+}
+
static int do_cpt_df_hist(XDR *xd, int fflags, int nlambda, df_history_t **dfhistPtr, FILE *list)
{
int ret = 0;
}
}
+ PullHistory *pullHist = observablesHistory->pullHistory.get();
+ int flagsPullHistory = 0;
+ if (pullHist != nullptr && (pullHist->numValuesInXSum > 0 || pullHist->numValuesInFSum > 0))
+ {
+ flagsPullHistory |= (1<<epullhPULL_NUMCOORDINATES);
+ flagsPullHistory |= ((1<<epullhPULL_NUMGROUPS) | (1<<epullhPULL_NUMVALUESINXSUM) | (1<<epullhPULL_NUMVALUESINFSUM));
+ }
+
int flags_dfh;
if (bExpanded)
{
eIntegrator, simulation_part, step, t, nppnodes,
{0}, npmenodes,
state->natoms, state->ngtc, state->nnhpres,
- state->nhchainlength, nlambda, state->flags, flags_eks, flags_enh, flags_dfh, flags_awhh,
+ state->nhchainlength, nlambda, state->flags, flags_eks, flags_enh,
+ flagsPullHistory, flags_dfh, flags_awhh,
nED, eSwapCoords
};
std::strcpy(headerContents.version, gmx_version());
if ((do_cpt_state(gmx_fio_getxdr(fp), state->flags, state, nullptr) < 0) ||
(do_cpt_ekinstate(gmx_fio_getxdr(fp), flags_eks, &state->ekinstate, nullptr) < 0) ||
(do_cpt_enerhist(gmx_fio_getxdr(fp), FALSE, flags_enh, enerhist, nullptr) < 0) ||
+ (doCptPullHist(gmx_fio_getxdr(fp), FALSE, flagsPullHistory, pullHist, StatePart::pullHistory, nullptr) < 0) ||
(do_cpt_df_hist(gmx_fio_getxdr(fp), flags_dfh, nlambda, &state->dfhist, nullptr) < 0) ||
(do_cpt_EDstate(gmx_fio_getxdr(fp), FALSE, nED, edsamhist, nullptr) < 0) ||
(do_cpt_awh(gmx_fio_getxdr(fp), FALSE, flags_awhh, state->awhHistory.get(), nullptr) < 0) ||
cp_error();
}
+ if (headerContents->flagsPullHistory)
+ {
+ if (observablesHistory->pullHistory == nullptr)
+ {
+ observablesHistory->pullHistory = gmx::compat::make_unique<PullHistory>();
+ }
+ ret = doCptPullHist(gmx_fio_getxdr(fp), TRUE,
+ headerContents->flagsPullHistory, observablesHistory->pullHistory.get(), StatePart::pullHistory, nullptr);
+ if (ret)
+ {
+ cp_error();
+ }
+ }
+
if (headerContents->file_version < 6)
{
gmx_fatal(FARGS, "Continuing from checkpoint files written before GROMACS 4.5 is not supported");
{
cp_error();
}
+ PullHistory pullHist = {};
+ ret = doCptPullHist(gmx_fio_getxdr(fp), TRUE,
+ headerContents.flagsPullHistory, &pullHist, StatePart::pullHistory, nullptr);
+ if (ret)
+ {
+ cp_error();
+ }
+
ret = do_cpt_df_hist(gmx_fio_getxdr(fp), headerContents.flags_dfh, headerContents.nlambda, &state->dfhist, nullptr);
if (ret)
{
ret = do_cpt_enerhist(gmx_fio_getxdr(fp), TRUE,
headerContents.flags_enh, &enerhist, out);
+ if (ret == 0)
+ {
+ PullHistory pullHist = {};
+ ret = doCptPullHist(gmx_fio_getxdr(fp), TRUE,
+ headerContents.flagsPullHistory, &pullHist, StatePart::pullHistory, out);
+ }
+
if (ret == 0)
{
ret = do_cpt_df_hist(gmx_fio_getxdr(fp),
tpxv_RemoveImplicitSolvation, /**< removed support for implicit solvation */
tpxv_PullPrevStepCOMAsReference, /**< Enabled using the COM of the pull group of the last frame as reference for PBC */
tpxv_MimicQMMM, /**< Inroduced support for MiMiC QM/MM interface */
+ tpxv_PullAverage, /**< Added possibility to output average pull force and position */
tpxv_Count /**< the total number of tpxv versions */
};
bRead, file_version, ePullOld, eGeomOld, dimOld);
}
}
+ if (file_version >= tpxv_PullAverage)
+ {
+ gmx_bool v;
+
+ v = pull->bXOutAverage;
+ gmx_fio_do_gmx_bool(fio, v);
+ pull->bXOutAverage = v;
+ v = pull->bFOutAverage;
+ gmx_fio_do_gmx_bool(fio, v);
+ pull->bFOutAverage = v;
+ }
}
pull->nstxout = get_eint(inp, "pull-nstxout", 50, wi);
pull->nstfout = get_eint(inp, "pull-nstfout", 50, wi);
pull->bSetPbcRefToPrevStepCOM = (get_eeenum(inp, "pull-pbc-ref-prev-step-com", yesno_names, wi) != 0);
+ pull->bXOutAverage = (get_eeenum(inp, "pull-xout-average", yesno_names, wi) != 0);
+ pull->bFOutAverage = (get_eeenum(inp, "pull-fout-average", yesno_names, wi) != 0);
printStringNoNewline(inp, "Number of pull groups");
pull->ngroup = get_eint(inp, "pull-ngroups", 1, wi);
printStringNoNewline(inp, "Number of pull coordinates");
#include "gromacs/mdtypes/md_enums.h"
#include "gromacs/mdtypes/mdatom.h"
#include "gromacs/mdtypes/observableshistory.h"
+#include "gromacs/mdtypes/pullhistory.h"
#include "gromacs/mdtypes/state.h"
#include "gromacs/pbcutil/mshift.h"
#include "gromacs/pbcutil/pbc.h"
{
restore_energyhistory_from_state(mdebin, observablesHistory->energyHistory.get());
}
- else if (observablesHistory->energyHistory != nullptr)
+ else
{
- /* We might have read an energy history from checkpoint.
- * As we are not appending, we want to restart the statistics.
- * Free the allocated memory and reset the counts.
+ if (observablesHistory->energyHistory != nullptr)
+ {
+ /* We might have read an energy history from checkpoint.
+ * As we are not appending, we want to restart the statistics.
+ * Free the allocated memory and reset the counts.
+ */
+ observablesHistory->energyHistory = {};
+ }
+ /* We might have read a pull history from checkpoint.
+ * We will still want to keep the statistics, so that the files
+ * can be joined and still be meaningful.
+ * This means that observablesHistory->pullHistory
+ * should not be reset.
*/
- observablesHistory->energyHistory = {};
}
if (ir->pull && ir->pull->bSetPbcRefToPrevStepCOM)
{
{
observablesHistory->energyHistory = compat::make_unique<energyhistory_t>();
}
+ if (observablesHistory->pullHistory == nullptr)
+ {
+ observablesHistory->pullHistory = compat::make_unique<PullHistory>();
+ }
/* Set the initial energy history in state by updating once */
update_energyhistory(observablesHistory->energyHistory.get(), mdebin);
}
inputrec->pull_work =
init_pull(fplog, inputrec->pull, inputrec,
&mtop, cr, &atomSets, inputrec->fepvals->init_lambda);
+ if (inputrec->pull->bXOutAverage || inputrec->pull->bFOutAverage)
+ {
+ initPullHistory(inputrec->pull_work, &observablesHistory);
+ }
if (EI_DYNAMICS(inputrec->eI) && MASTER(cr))
{
init_pull_output_files(inputrec->pull_work,
PI("pull-nstxout", pull->nstxout);
PI("pull-nstfout", pull->nstfout);
PS("pull-pbc-ref-prev-step-com", EBOOL(pull->bSetPbcRefToPrevStepCOM));
+ PS("pull-xout-average", EBOOL(pull->bXOutAverage));
+ PS("pull-fout-average", EBOOL(pull->bFOutAverage));
PI("pull-ngroups", pull->ngroup);
for (g = 0; g < pull->ngroup; g++)
{
/*
* This file is part of the GROMACS molecular simulation package.
*
- * Copyright (c) 2017, by the GROMACS development team, led by
+ * Copyright (c) 2017,2018, 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/mdtypes/edsamhistory.h"
#include "gromacs/mdtypes/energyhistory.h"
+#include "gromacs/mdtypes/pullhistory.h"
#include "gromacs/mdtypes/swaphistory.h"
ObservablesHistory::ObservablesHistory() = default;
#include <memory>
class energyhistory_t;
+class PullHistory;
struct edsamhistory_t;
struct swaphistory_t;
//! History for energy observables, used for output only
std::unique_ptr<energyhistory_t> energyHistory;
+ //! History for pulling observables, used for output only
+ std::unique_ptr<PullHistory> pullHistory;
+
//! Essential dynamics and flooding history
std::unique_ptr<edsamhistory_t> edsamHistory;
gmx_bool bSetPbcRefToPrevStepCOM; /**< Use the COM of each group from the previous step as reference */
int nstxout; /**< Output interval for pull x */
int nstfout; /**< Output interval for pull f */
+ bool bXOutAverage; /**< Write the average coordinate during the output interval */
+ bool bFOutAverage; /**< Write the average force during the output interval */
t_pull_group *group; /**< groups to pull/restrain/etc/ */
t_pull_coord *coord; /**< the pull coordinates */
--- /dev/null
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2016,2017,2018, 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.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+
+/*! \libinternal \file
+ *
+ *
+ * \brief
+ * This file contains datatypes for pull statistics history.
+ *
+ * \author Magnus Lundborg, Berk Hess
+ *
+ * \inlibraryapi
+ * \ingroup module_mdtypes
+ */
+
+#ifndef GMX_MDLIB_PULLHISTORY_H
+#define GMX_MDLIB_PULLHISTORY_H
+
+#include <vector>
+
+//! \cond INTERNAL
+
+//! \brief Contains the sum of coordinate observables to enable calculation of the average of pull data.
+class PullCoordinateHistory
+{
+ public:
+ double value; //!< The sum of the current value of the coordinate, units of nm or rad.
+ double valueRef; //!< The sum of the reference value of the coordinate, units of nm or rad.
+ double scalarForce; //!< The sum of the scalar force of the coordinate.
+ dvec dr01; //!< The sum of the direction vector of group 1 relative to group 0.
+ dvec dr23; //!< The sum of the direction vector of group 3 relative to group 2.
+ dvec dr45; //!< The sum of the direction vector of group 5 relative to group 4.
+ dvec dynaX; //!< The sum of the coordinate of the dynamic groups for geom=cylinder.
+
+ //! Constructor
+ PullCoordinateHistory() : value (0),
+ valueRef (0),
+ scalarForce(0),
+ dr01 (),
+ dr23 (),
+ dr45 (),
+ dynaX ()
+ {
+ }
+};
+
+//! \brief Contains the sum of group observables to enable calculation of the average of pull data.
+class PullGroupHistory
+{
+ public:
+ dvec x; //!< The sum of the coordinates of the group.
+
+ //! Constructor
+ PullGroupHistory() : x ()
+ {
+ }
+};
+
+
+//! \brief Pull statistics history, to allow output of average pull data.
+class PullHistory
+{
+ public:
+ int numValuesInXSum; //!< Number of values of the coordinate values in the pull sums.
+ int numValuesInFSum; //!< Number of values in the pull force sums.
+ std::vector<PullCoordinateHistory> pullCoordinateSums; //!< The container of the sums of the values of the pull coordinate, also contains the scalar force.
+ std::vector<PullGroupHistory> pullGroupSums; //!< The container of the sums of the values of the pull group.
+
+ //! Constructor
+ PullHistory() : numValuesInXSum(0),
+ numValuesInFSum(0)
+ {
+ }
+};
+
+//! \endcond
+
+#endif
#include <cstdio>
#include "gromacs/commandline/filenm.h"
+#include "gromacs/compat/make_unique.h"
#include "gromacs/fileio/gmxfio.h"
#include "gromacs/fileio/xvgr.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/mdtypes/observableshistory.h"
+#include "gromacs/mdtypes/pullhistory.h"
#include "gromacs/pulling/pull.h"
#include "gromacs/utility/cstringutil.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
#include "pull_internal.h"
}
}
-static void pull_print_group_x(FILE *out, const ivec dim,
- const pull_group_work_t *pgrp)
+static void addToPullxHistory(pull_t *pull)
{
- int m;
+ GMX_RELEASE_ASSERT(pull->coordForceHistory, "Pull history does not exist.");
+ pull->coordForceHistory->numValuesInXSum++;
+ for (size_t c = 0; c < pull->coord.size(); c++)
+ {
+ const pull_coord_work_t &pcrd = pull->coord[c];
+ PullCoordinateHistory &pcrdHistory = pull->coordForceHistory->pullCoordinateSums[c];
- for (m = 0; m < DIM; m++)
+ pcrdHistory.value += pcrd.spatialData.value;
+ pcrdHistory.valueRef += pcrd.value_ref;
+
+ for (int m = 0; m < DIM; m++)
+ {
+ pcrdHistory.dr01[m] += pcrd.spatialData.dr01[m];
+ pcrdHistory.dr23[m] += pcrd.spatialData.dr23[m];
+ pcrdHistory.dr45[m] += pcrd.spatialData.dr45[m];
+ }
+ if (pcrd.params.eGeom == epullgCYL)
+ {
+ for (int m = 0; m < DIM; m++)
+ {
+ pcrdHistory.dynaX[m] += pull->dyna[c].x[m];
+ }
+ }
+ }
+ for (size_t g = 0; g < pull->group.size(); g++)
{
- if (dim[m])
+ PullGroupHistory &pgrpHistory = pull->coordForceHistory->pullGroupSums[g];
+ for (int m = 0; m < DIM; m++)
+ {
+ pgrpHistory.x[m] += pull->group[g].x[m];
+ }
+ }
+}
+
+static void addToPullfHistory(pull_t *pull)
+{
+ GMX_RELEASE_ASSERT(pull->coordForceHistory, "Pull history does not exist.");
+ pull->coordForceHistory->numValuesInFSum++;
+ for (size_t c = 0; c < pull->coord.size(); c++)
+ {
+ const pull_coord_work_t &pcrd = pull->coord[c];;
+ PullCoordinateHistory &pcrdHistory = pull->coordForceHistory->pullCoordinateSums[c];
+
+ pcrdHistory.scalarForce += pcrd.scalarForce;
+ }
+}
+
+static void pullResetHistory(PullHistory *history,
+ bool resetXHistory,
+ bool resetFHistory)
+{
+ if (resetXHistory)
+ {
+ history->numValuesInXSum = 0;
+
+ for (PullCoordinateHistory &pcrdHistory : history->pullCoordinateSums)
+ {
+ pcrdHistory.value = 0;
+ pcrdHistory.valueRef = 0;
+
+ clear_dvec(pcrdHistory.dr01);
+ clear_dvec(pcrdHistory.dr23);
+ clear_dvec(pcrdHistory.dr45);
+ clear_dvec(pcrdHistory.dynaX);
+ }
+
+ for (PullGroupHistory &pgrpHistory : history->pullGroupSums)
+ {
+ clear_dvec(pgrpHistory.x);
+ }
+ }
+ if (resetFHistory)
+ {
+ history->numValuesInFSum = 0;
+ for (PullCoordinateHistory &pcrdHistory : history->pullCoordinateSums)
{
- fprintf(out, "\t%g", pgrp->x[m]);
+ pcrdHistory.scalarForce = 0;
}
}
}
-static void pull_print_coord_dr_components(FILE *out, const ivec dim, const dvec dr)
+static void pull_print_coord_dr_components(FILE *out, const ivec dim, const dvec dr,
+ const int numValuesInSum)
{
for (int m = 0; m < DIM; m++)
{
if (dim[m])
{
- fprintf(out, "\t%g", dr[m]);
+ fprintf(out, "\t%g", dr[m] / numValuesInSum);
}
}
}
-static void pull_print_coord_dr(FILE *out, const pull_coord_work_t *pcrd,
- gmx_bool bPrintRefValue,
- gmx_bool bPrintComponents)
+template <typename T>
+static void pull_print_coord_dr(FILE *out,
+ const pull_params_t &pullParams,
+ const t_pull_coord &coordParams,
+ const T &pcrdData,
+ double referenceValue,
+ const int numValuesInSum)
{
- double unit_factor = pull_conversion_factor_internal2userinput(&pcrd->params);
+ const double unit_factor = pull_conversion_factor_internal2userinput(&coordParams);
- fprintf(out, "\t%g", pcrd->spatialData.value*unit_factor);
+ fprintf(out, "\t%g", pcrdData.value*unit_factor/numValuesInSum);
- if (bPrintRefValue)
+ if (pullParams.bPrintRefValue && coordParams.eType != epullEXTERNAL)
{
- fprintf(out, "\t%g", pcrd->value_ref*unit_factor);
+ fprintf(out, "\t%g", referenceValue*unit_factor/numValuesInSum);
}
- if (bPrintComponents)
+ if (pullParams.bPrintComp)
{
- pull_print_coord_dr_components(out, pcrd->params.dim, pcrd->spatialData.dr01);
- if (pcrd->params.ngroup >= 4)
+ pull_print_coord_dr_components(out, coordParams.dim, pcrdData.dr01, numValuesInSum);
+ if (coordParams.ngroup >= 4)
{
- pull_print_coord_dr_components(out, pcrd->params.dim, pcrd->spatialData.dr23);
+ pull_print_coord_dr_components(out, coordParams.dim, pcrdData.dr23, numValuesInSum);
}
- if (pcrd->params.ngroup >= 6)
+ if (coordParams.ngroup >= 6)
{
- pull_print_coord_dr_components(out, pcrd->params.dim, pcrd->spatialData.dr45);
+ pull_print_coord_dr_components(out, coordParams.dim, pcrdData.dr45, numValuesInSum);
}
}
}
-static void pull_print_x(FILE *out, struct pull_t *pull, double t)
+static void pull_print_x(FILE *out, pull_t *pull, double t)
{
fprintf(out, "%.4f", t);
for (size_t c = 0; c < pull->coord.size(); c++)
{
- pull_coord_work_t *pcrd;
+ const pull_coord_work_t &pcrd = pull->coord[c];
+ int numValuesInSum = 1;
+ const PullCoordinateHistory *pcrdHistory = nullptr;
- pcrd = &pull->coord[c];
+ if (pull->bXOutAverage)
+ {
+ pcrdHistory = &pull->coordForceHistory->pullCoordinateSums[c];
- pull_print_coord_dr(out, pcrd,
- pull->params.bPrintRefValue && pcrd->params.eType != epullEXTERNAL,
- pull->params.bPrintComp);
+ numValuesInSum = pull->coordForceHistory->numValuesInXSum;
+ pull_print_coord_dr(out, pull->params, pcrd.params,
+ *pcrdHistory, pcrdHistory->valueRef,
+ numValuesInSum);
+ }
+ else
+ {
+ pull_print_coord_dr(out, pull->params, pcrd.params,
+ pcrd.spatialData, pcrd.value_ref,
+ numValuesInSum);
+ }
if (pull->params.bPrintCOM)
{
- if (pcrd->params.eGeom == epullgCYL)
+ if (pcrd.params.eGeom == epullgCYL)
{
- pull_print_group_x(out, pcrd->params.dim, &pull->dyna[c]);
+ for (int m = 0; m < DIM; m++)
+ {
+ if (pcrd.params.dim[m])
+ {
+ /* This equates to if (pull->bXOutAverage) */
+ if (pcrdHistory)
+ {
+ fprintf(out, "\t%g", pcrdHistory->dynaX[m] / numValuesInSum);
+ }
+ else
+ {
+ fprintf(out, "\t%g", pull->dyna[c].x[m]);
+ }
+ }
+ }
}
else
{
- pull_print_group_x(out, pcrd->params.dim,
- &pull->group[pcrd->params.group[0]]);
+ for (int m = 0; m < DIM; m++)
+ {
+ if (pcrd.params.dim[m])
+ {
+ if (pull->bXOutAverage)
+ {
+ fprintf(out, "\t%g", pull->coordForceHistory->pullGroupSums[pcrd.params.group[0]].x[m] / numValuesInSum);
+ }
+ else
+ {
+ fprintf(out, "\t%g", pull->group[pcrd.params.group[0]].x[m]);
+ }
+ }
+ }
}
- for (int g = 1; g < pcrd->params.ngroup; g++)
+ for (int g = 1; g < pcrd.params.ngroup; g++)
{
- pull_print_group_x(out, pcrd->params.dim, &pull->group[pcrd->params.group[g]]);
+ for (int m = 0; m < DIM; m++)
+ {
+ if (pcrd.params.dim[m])
+ {
+ if (pull->bXOutAverage)
+ {
+ fprintf(out, "\t%g", pull->coordForceHistory->pullGroupSums[pcrd.params.group[g]].x[m] / numValuesInSum);
+ }
+ else
+ {
+ fprintf(out, "\t%g", pull->group[pcrd.params.group[g]].x[m]);
+ }
+ }
+ }
}
}
}
fprintf(out, "\n");
+
+ if (pull->bXOutAverage)
+ {
+ pullResetHistory(pull->coordForceHistory, TRUE, FALSE);
+ }
}
static void pull_print_f(FILE *out, const pull_t *pull, double t)
{
fprintf(out, "%.4f", t);
- for (const pull_coord_work_t &coord : pull->coord)
+ if (pull->bFOutAverage)
{
- fprintf(out, "\t%g", coord.scalarForce);
+ for (size_t c = 0; c < pull->coord.size(); c++)
+ {
+ fprintf(out, "\t%g", pull->coordForceHistory->pullCoordinateSums[c].scalarForce / pull->coordForceHistory->numValuesInFSum);
+ }
+ }
+ else
+ {
+ for (const pull_coord_work_t &coord : pull->coord)
+ {
+ fprintf(out, "\t%g", coord.scalarForce);
+ }
}
fprintf(out, "\n");
+
+ if (pull->bFOutAverage)
+ {
+ pullResetHistory(pull->coordForceHistory, FALSE, TRUE);
+ }
}
void pull_print_output(struct pull_t *pull, int64_t step, double time)
{
GMX_ASSERT(pull->numExternalPotentialsStillToBeAppliedThisStep == 0, "pull_print_output called before all external pull potentials have been applied");
- if ((pull->params.nstxout != 0) && (step % pull->params.nstxout == 0))
+ if (pull->params.nstxout != 0)
{
- pull_print_x(pull->out_x, pull, time);
+ /* Do not update the average if the number of observations already equal (or are
+ * higher than) what should be in each average output. This can happen when
+ * appending to a file from a checkpoint, which would otherwise include the
+ * last value twice.*/
+ if (pull->bXOutAverage && !pull->coord.empty() && pull->coordForceHistory->numValuesInXSum < pull->params.nstxout)
+ {
+ addToPullxHistory(pull);
+ }
+ if (step % pull->params.nstxout == 0)
+ {
+ pull_print_x(pull->out_x, pull, time);
+ }
}
- if ((pull->params.nstfout != 0) && (step % pull->params.nstfout == 0))
+ if (pull->params.nstfout != 0)
{
- pull_print_f(pull->out_f, pull, time);
+ /* Do not update the average if the number of observations already equal (or are
+ * higher than) what should be in each average output. This can happen when
+ * appending to a file from a checkpoint, which would otherwise include the
+ * last value twice.*/
+ if (pull->bFOutAverage && !pull->coord.empty() && pull->coordForceHistory->numValuesInFSum < pull->params.nstfout)
+ {
+ addToPullfHistory(pull);
+ }
+ if (step % pull->params.nstfout == 0)
+ {
+ pull_print_f(pull->out_f, pull, time);
+ }
}
}
if (bCoord)
{
sprintf(buf, "Position (nm%s)", pull->bAngle ? ", deg" : "");
- xvgr_header(fp, "Pull COM", "Time (ps)", buf,
- exvggtXNY, oenv);
+ if (pull->bXOutAverage)
+ {
+ xvgr_header(fp, "Pull Average COM", "Time (ps)", buf,
+ exvggtXNY, oenv);
+ }
+ else
+ {
+ xvgr_header(fp, "Pull COM", "Time (ps)", buf,
+ exvggtXNY, oenv);
+ }
}
else
{
sprintf(buf, "Force (kJ/mol/nm%s)", pull->bAngle ? ", kJ/mol/rad" : "");
- xvgr_header(fp, "Pull force", "Time (ps)", buf,
- exvggtXNY, oenv);
+ if (pull->bFOutAverage)
+ {
+ xvgr_header(fp, "Pull Average force", "Time (ps)", buf,
+ exvggtXNY, oenv);
+ }
+ else
+ {
+ xvgr_header(fp, "Pull force", "Time (ps)", buf,
+ exvggtXNY, oenv);
+ }
}
/* With default mdp options only the actual coordinate value is printed (1),
FALSE, continuationOptions);
}
}
+
+void initPullHistory(pull_t *pull,
+ ObservablesHistory *observablesHistory)
+{
+ GMX_RELEASE_ASSERT(pull, "Need a valid pull object");
+
+ if (observablesHistory == nullptr)
+ {
+ pull->coordForceHistory = nullptr;
+ return;
+ }
+ /* If pull->coordForceHistory is already set we are starting from a checkpoint. Do not reset it. */
+ if (observablesHistory->pullHistory == nullptr)
+ {
+ observablesHistory->pullHistory = gmx::compat::make_unique<PullHistory>();
+ pull->coordForceHistory = observablesHistory->pullHistory.get();
+ pull->coordForceHistory->numValuesInXSum = 0;
+ pull->coordForceHistory->numValuesInFSum = 0;
+ pull->coordForceHistory->pullCoordinateSums.resize(pull->coord.size());
+ pull->coordForceHistory->pullGroupSums.resize(pull->group.size());
+ }
+ else
+ {
+ pull->coordForceHistory = observablesHistory->pullHistory.get();
+ }
+}
struct pull_t;
struct ContinuationOptions;
struct gmx_output_env_t;
+struct ObservablesHistory;
struct t_filenm;
/*! \brief Set up and open the pull output files, when requested.
*/
void pull_print_output(pull_t *pull, int64_t step, double time);
+/*! \brief Allocate and initialize pull work history (for average pull output) and set it in a pull work struct
+ *
+ * \param pull The pull work struct
+ * \param observablesHistory Container of history data, e.g., pull history.
+ */
+void initPullHistory(pull_t *pull, ObservablesHistory *observablesHistory);
+
#endif
const bool computeVirial = (force->computeVirial_ && MASTER(cr));
for (size_t c = 0; c < pull->coord.size(); c++)
{
+ pull_coord_work_t *pcrd;
+ pcrd = &pull->coord[c];
+
/* For external potential the force is assumed to be given by an external module by a call to
apply_pull_coord_external_force */
- if (pull->coord[c].params.eType == epullCONSTRAINT || pull->coord[c].params.eType == epullEXTERNAL)
+ if (pcrd->params.eType == epullCONSTRAINT || pcrd->params.eType == epullEXTERNAL)
{
continue;
}
}
}
- pull->bPotential = FALSE;
- pull->bConstraint = FALSE;
- pull->bCylinder = FALSE;
- pull->bAngle = FALSE;
+ pull->bPotential = FALSE;
+ pull->bConstraint = FALSE;
+ pull->bCylinder = FALSE;
+ pull->bAngle = FALSE;
+ pull->bXOutAverage = pull_params->bXOutAverage;
+ pull->bFOutAverage = pull_params->bFOutAverage;
GMX_RELEASE_ASSERT(pull->group[0].params.nat == 0, "pull group 0 is an absolute reference group and should not contain atoms");
pull->group[0].x_prev_step[XX] = NAN;
static const int c_pullMaxNumLocalAtomsSingleThreaded = 1;
#endif
+class PullHistory;
+
enum {
epgrppbcNONE, epgrppbcREFAT, epgrppbcCOS, epgrppbcPREVSTEPCOM
};
std::vector<pull_coord_work_t> coord; /* The pull group param and work data */
/* Global dynamic data */
- gmx_bool bSetPBCatoms; /* Do we need to set x_pbc for the groups? */
+ gmx_bool bSetPBCatoms; /* Do we need to set x_pbc for the groups? */
+
+ int nthreads; /* Number of threads used by the pull code */
+ std::vector<ComSums> comSums; /* Work array for summing for COM, 1 entry per thread */
+
+ pull_comm_t comm; /* Communication parameters, communicator and buffers */
- int nthreads; /* Number of threads used by the pull code */
- std::vector<ComSums> comSums; /* Work array for summing for COM, 1 entry per thread */
+ FILE *out_x; /* Output file for pull data */
+ FILE *out_f; /* Output file for pull data */
- pull_comm_t comm; /* Communication parameters, communicator and buffers */
+ bool bXOutAverage; /* Output average pull coordinates */
+ bool bFOutAverage; /* Output average pull forces */
- FILE *out_x; /* Output file for pull data */
- FILE *out_f; /* Output file for pull data */
+ PullHistory *coordForceHistory; /* Pull coordinate and force history */
/* The number of coordinates using an external potential */
int numCoordinatesWithExternalPotential;
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff